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(ABSTRACT) 

Design for prevention of aeroelastic instability (that is, the critical speeds 
leading to aeroelastic instability lie outside the operating range) is an integral part 
of the wing design process. Availability of the sensitivity derivatives of the various 
critical speeds with respect to shape parameters of the wing could be very useful 
to a designer in the initial design phase, when several design changes are made and 
the shape of the final configuration is not yet frozen. These derivatives are also 
indispensable for a gradient-based optimization with aeroelastic constraints. 

In this study, flutter characteristic of a typical section in subsonic com- 
pressible flow is examined using a state-space unsteady aerodynamic represen- 
tation. The sensitivity of the flutter speed of the typical section with respect to 
its mass and stiffness parameters, namely, mass ratio, static unbalance, radius 
of gyration, bending frequency and torsional frequency is calculated analytically. 
A strip-theory formulation is newly developed to represent the unsteady aerody- 
namic forces on a wing. This is coupled with an equivalent plate structural model 
based on a Rayleigh- Ritz formulation and the aeroelastic equations .are solved as 
an eigenvalue problem to determine the critical speed of the wing. The sensitivity 



of divergence and flutter speeds to shape parameters, namely, aspect ratio, area, 
taper ratio and sweep angle are computed analytically. The aeroelastic equations 
are also integrated with respect to time using the Wilson-0 method at different 
values of freestream speed, to observe the aeroelastic phenomena in real time. 

Flutter analysis of the wing is also carried out using a lifting-surface subsonic 
kernel function aerodynamic theory (FAST) and an equivalent plate structural 
model. The flutter speed is obtained using a V-g type of solution. The novel 
method of automatic differentiation using ADIFOR is implemented to generate 
exact derivatives of the flutter speed with respect to shape and modal parameters 
of the wing. 

Finite element modeling of the wing is done using NASTRAN so that wing 
structures made of spars and ribs and top and bottom wing skins could be analyzed. 
The free vibration modes of the wing obtained from NASTRAN axe input into 
FAST to compute the flutter speed. The derivatives of flutter speed with respect 
to shape parameters are computed using a combination of central difference scheme 
and ADIFOR and the sensitivity to modal parameters is calculated using ADIFOR. 

An equivalent plate model which incorporates first-order shear deformation 
theory is then examined so it can be used to model thick wings, where shear 
deformations are important. The sensitivity of natural frequencies to changes 
in shape parameters is obtained using ADIFOR. A simple optimization effort is 
made towards obtaining a minimum weight design of the wing, subject to flutter 
constraints, lift requirement constraints for level flight and side constraints on the 
planform parameters of the wing using the IMSL subroutine NCONG, which uses 
successive quadratic programming. 
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CHAPTER 1 
INTRODUCTION 


1.1 BACKGROUND 

Since aircraft structures are flexible and they deform under the aerodynamic 
loads encountered during flight, studies on structure- fluid interactions and the re- 
sponse of the aircraft structure under these conditions are very important. The 
interaction between the elastic forces and the aerodynamic forces on a flexible 
structure in motion alters its vibration characteristics. In an aircraft structure, for 
example, the lift forces and moments on the wing are modified by the elastic defor- 
mation of the structure and the new airloads, in turn, produce a new deformation 
pattern. The aerodynamic forces causing the deformation are proportional to the 
square of the velocity, whereas the elastic restoring forces are proportional to the 
stiffness of the structure. Therefore, at speeds above the critical speeds, when 
the aerodynamic forces exceed the elastic restoring forces or when the structure 
undergoes oscillations of increasing amplitude forced by the oscillatory unsteady 
airloads, aeroelastic instabilities [1—3] like divergence and flutter occur. 

In a multidisciplinary design environment, the problems involving interac- 
tions between two or more disciplines have to be addressed. Design for aeroelastic 
stability of an aircraft is an integral part of the wing design process. It is impera- 
tive during aircraft design that the critical speeds leading to instability lie outside 
the operating range of the aircraft. The structural and aerodynamic character- 
istics of a wing are functions of its shape and hence the aeroelastic response is 
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sensitive to changes in shape parameters. Sensitivity analysis is an important tool 
which yields information about the dependence of the aeroelastic instability on the 
design parameters of the wing. 

Sensitivity derivatives of the aeroelastic response of a wing with respect to 
its shape parameters are very useful for preliminary design purposes. In concep- 
tual and preliminary design, shape variations of the airplane should be considered, 
before the shape of the configuration is frozen. The generalized unsteady aerody- 
namic forces need to be updated with changes in wing planform and have to be 
evaluated at several values of reduced frequencies before an aeroelastic analysis can 
be performed. However, if the sensitivity derivatives are computed at the baseline 
configuration, it gives a linear approximation to the aeroelastic response curve and 
can be used for small changes in the shape design parameters from the baseline, 
without having to perform a reanalysis at the new configuration. 

Sensitivity derivatives are also of great importance in multidisciplinary design 
optimization. It is very beneficial in regards to time and cost, when an integrated 
multidisciplinary design and synthesis approach is followed, where all the relevant 
disciplines are considered simultaneously [4—6] . The systems that are decompos- 
able into top-down hierarchy of engineering disciplines and subsystems may be 
optimized by multilevel procedure made up of suboptimizations performed concur- 
rently at each level of the hierarchy and linked by optimum sensitivity derivative 
information [7]. In a gradient-based optimization with aeroelastic constraints, the 
sensitivity derivatives of the aeroelastic response of a wing will comprise the local 
sensitivity information at the subsystem level. 
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The sensitivity derivatives of a system can be found using either analytical 
or finite difference methods. Analytical sensitivity analysis has found increased 
interest in engineering design as it eliminates uncertainity in the choice of step size 
needed in the finite difference method. The step size if too large leads to truncation 
errors and if too small leads to ill-conditioning. 

1.2 LITERATURE REVIEW 

1.2.1 AEROELASTIC ANALYSIS 

The various methods that are used for aeroelastic analysis of a typical section 
or a wing differ in the prediction of aerodynamic loads. For many years, aeroelastic 
analysis has been performed using linearized aerodynamic and structural theories 
[8-12]. To generate sensitivity derivatives for use in preliminary design or an opti- 
mization routine, reasonably quick and accurate methods for aeroelastic analysis 
are desired. The lift and moment predictions on an airfoil undergoing harmonic 
motion have been obtained by Theodorsen [13]. Several CFD methods which are 
used to determine the transonic flowfield around two dimensional airfoils are listed 
by Ballhaus and Bridgeman [14] and used in [15-17]. 

In recent years, considerable effort has been made to integrate the aerody- 
namic, structural and control aspects of the design of an aircraft. Livne et al 
[18] comment that for integrated multidisciplinary wing synthesis, where design 
for aeroservoelastic stability is an objective, it is required to represent the aeroe- 
lastic equations of motion in Linear Time Invariant (LTI) state-space form. The 
unsteady aerodynamic loads on the wing can be represented in a state-space form, 
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thereby adding only a small number of states to the mathematical model of the 
aeroservoelastic system. 

The unsteady aerodynamic loads on a typical section in response to an arbi- 
trary forcing can be represented through the use of indicial (step) response func- 
tions (19). The indicial response method is a fundamental approach to the problem, 
and affords considerable insight into the physical aspects of unsteady airfoil flow. 
One main advantage of the approach is that when the indicial response to a par- 
ticular forcing mode is known, e.g., that due to angle of attack or pitch rate, the 
cumulative response to an arbitrary forcing can be obtained in the time-domain 
by means of Duhamel superposition. Indicial response for an incompressible flow 
was obtained theoretically by Wagner [20], Jones [21] used a two-pole exponential 
approximation to the Wagner function. Venkatesan and Friedmann [22] have given 
a three-pole indicial response function that can express the Theodorsen’s function 
over the entire reduced frequency range. 

Leishman and Nguyen [23] have represented the aerodynamic indicial re- 
sponse functions for compressible flow by upto three-pole approximations, the 
response consisting of two parts, one due to non-circulatory loading and the other 
due to circulatory loading. This has advantages over the CFD-based methods 
in the sense that the CFD methods are in general computationally very expen- 
sive, The state equations describing the unsteady aerodynamic response can be 
obtained by direct application of Laplace transforms to these indicial response 
functions [24,25]. State-space representation of aerodynamic characteristics of an 
aircraft at high angles of attack considering flow separation and flow with vortex 
breakdown has been presented by Goman and Khrabrov [26]. 
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Aerodynamic methods based on strip-theory formulations [27,28] and ad- 
vanced codes such as XTRAN3S [29] and CAP-TSD [30], which use the transonic 
small disturbance equation are currently being used for aeroelastic analysis. Lin- 
ear lifting-surface theories such as the Doublet Lattice method [31,32] and Kernel 
function methods [33,34] give adequate prediction of aerodynamic loads required 
for aeroelastic analysis in the subsonic and supersonic regimes. 

1.2.2 SENSITIVITY ANALYSIS 

Sensitivity analysis is becoming an important design tool in engineering de- 
sign applications. It was first recognized as a useful tool for assessing the effects 
of changing parameters in mathematical models of control systems. The gradient 
based mathematical programming method used in optimal control and structural 
optimization furthered the development of sensitivity derivatives, because sen- 
sitivity derivatives are used in search directions to find optimum solution [35]. 
Sensitivity analysis has also become a versatile design tool, rather than just an 
instrument of optimization programs [36]. Sobieski [37] discusses in detail about 
the System Design Derivatives which help in understanding the effect a particular 
design variable would have on the desired performance of the system, if it were 
perturbed by a small percentage from its original value. 

Adelman and Haftka [36] have shown that structural sensitivity analysis has 
been available for over two decades. Rudisill and Bhatia [38] developed expressions 
for the analytical derivatives of the eigenvalues, reduced frequency and flutter speed 
with respect to structural parameters for use in minimizing the total mass. Murthy 
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and Haftka [39] have presented a survey of methods for calculating sensitivity of 
general eigenproblems. Eigenvalue and eigenvector derivatives for general matrices 
are discussed and various approximation methods for eigenvalues are proposed and 
evaluated by Murthy [40], Pedersen and Seyranian [41] examined the change in 
flutter load as a function of change in stiffness, mass, boundary conditions or load 
distribution. They performed sensitivity analysis without resorting to any new 
eigenvalue analysis. The solution to the main and an adjoint problem provide all 
the necessary information for evaluating sensitivities. Their paper mainly focused 
on column and beam critical load distributions. 

Hawk and Bristow [42] developed aerodynamic sensitivity analysis capabili- 
ties in subcritical compressible flow. They first analyzed a baseline configuration, 
and then calculated a matrix containing partial derivatives of the potential at each 
control point with respect to each known geometric parameter by applying a first 
order expansion to the baseline configuration. The matrix of partial derivatives is 
used in each iteration cycle to analyze the perturbed geometry. However, this anal- 
ysis only handles chordwise perturbation distributions, such as changes in camber, 
thickness and twist. Another approach has been presented by Yates [43] that con- 
siders general geometric variations, including planform, and subsonic, sonic and 
supersonic unsteady, nonplanar lifting-surface theory. 

Recently, Livne et al [44] applied an equivalent plate structural modeling, 
which includes transverse shear, to an HSCT wing. Simple polynomials were used 
for Ritz functions and depth and thickness distributions. The derivatives of the 
stiffness and mass matrices were obtained analytically with respect to the shape 


6 



variables of the wing. Livne [45] observed that as higher order polynomials are 
used for better modeling of the structure, the more sensitive is the finite difference 
derivative to the step-size used and in some cases, it is impossible to obtain any 
valuable information by finite differences, 

Barthelemy and Bergen [46] explored the analytical shape sensitivity deriva- 
tives of the wing’s aeroelastic characteristics, such as section lift, angle of attack, 
rolling moment, induced drag and divergence dynamic pressure, for subsonic sub- 
critical flow, with respect to geometric parameters. Results showed the charac- 
teristics nonlinearity to be small enough to be well approximated by sensitivity 
based linear approximations. These approximations are valid within a range that 
is useful to designers in the initial design phase. 

Kapania [47] has obtained sensitivity derivatives of the flutter speed of a two 
dimensional airfoil in incompressible flow with respect to the mass and stiffness 
parameters. Kapania, Bergen and Barthelemy [28] have obtained the shape sen- 
sitivity derivatives of the flutter response of a laminated wing in incompressible 
flow. In this work, Yates’ modified strip analysis [48] was used for the aerodynamic 
model in conjunction with Giles’ equivalent plate analysis [49,50] for the structural 
model. 

1.2.3 MULTIDISCIPLINARY OPTIMIZATION 

Research. carried out on structural design optimization under various static 
and dynamic behavior constraints including flutter and static aeroelastic con- 
straints has been presented in [51-65]. Rudisill and Bhatia [38] used projected 
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gradient search to arrive at a relative optimum design of the structure for a speci- 
fied minimum flutter velocity. Further, they made use of second derivatives of the 
flutter velocity with respect to structural design parameters to reduce design cycles 
[66]. Haftka [67] has written a survey paper on application of structural optimiza- 
tion techniques to problems of design under aeroelastic constraints. Haftka and 
Yates [68] adopted a method of avoiding repetition of calculations which do not 
depend on the structural sizes, when repeated flutter calculations have to be made. 
Haftka and Prasad [69] resorted to computing vibration modes only periodically. 
For intermediate calculations, old modes were used as generalized coordinates to 
calculate approximate new modes, thus mostly avoiding the expensive direct vi- 
bration modes calculation. 

Sensitivity derivatives are of great importance in integrated multidisciplinary 
design optimization of aircrafts. More attention is directed towards multidisci- 
plinary structural/aerodynamic synthesis of wings [70-72]. Karpel [73] used a 
gradient-based constrained optimization on a composite active-flexible wing to 
achieve aircraft performance requirements and sufficient flutter and control sta- 
bility margins with a minimum weight penalty and without violating the design 
constraints. The sensitivity derivatives of the flutter dynamic pressure, control 
stability margins and control effectiveness with respect to structural and control 
design variables were obtained analytically. 

Hajela et al [74] applied Sobieski’s Global Sensitivity Equations (GSE) in an 
aircraft synthesis problem where the constraints involved the coupled disciplines 
of structures, aerodynamics and flight mechanics. The coupled system was repre- 
sented by smaller subsystems and the total behavior sensitivities were determined 
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by applying the GSE method. Barthelemy et al [75] discuss a multidisciplinary 
design optimization method applied to a supersonic transport wing. Aerodynamic 
and structural disciplines are integrated for a minimum weight design under static 
aeroelastic constraints. The authors point out that as the number of dependent 
variables in each discipline becomes large, the calculation of the finite difference 
derivatives contributes substantially to the total optimization cost. 

1.3 OBJECTIVES 

The primary objective of this research is to calculate sensitivity derivatives 
of aeroelastic responses at the critical condition (i.e., flutter and divergence) of a 
wing for use in preliminary design or a multidisciplinary optimization with aeroe- 
lastic constraints. Though advanced finite element structural models and com- 
putational fluid dynamics (CFD) models can be adopted for a detailed analysis, 
simpler structural and aerodynamic models are desired for aeroelastic analysis at 
the preliminary stages of design. Both equivalent plate and finite element models 
are used in this study for structural modeling and strip-theory and lifting-surface 
theory are used for the unsteady aerodynamics. Since analytic expressions exist for 
computing the derivatives of the eigenvalues of a general matrix, if the aeroelastic 
stability problem is solved as an eigenvalue problem, then sensitivities of the criti- 
cal speeds with respect to various wing parameters must be computed analytically 
whenever possible so that finite difference methods do not have to be resorted to. 
With this objective in mind, the sensitivity derivatives of the critical speeds have 
been calculated analytically as well as using automatic differentiation. 
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In Chapter 2, the onset of aeroelastic instability of a typical section in sub- 
sonic compressible flow is determined using a state-space unsteady aerodynamics 
representation. The sensitivity of flutter speed with respect to the mass and stiff- 
ness parameters are calculated analytically. A strip-theory aerodynamic formula- 
tion based on this state-space model is coupled with an equivalent plate structural 
model to determine the critical speed of a wing in Chapter 3. The sensitivity 
of divergence and flutter speeds to shape parameters, namely, aspect ratio, area, 
taper ratio and sweep angle are computed analytically. The aeroelastic equations 
are also integrated with respect to time using the Wilson-0 method. Flutter anal- 
ysis is carried out using a lifting-surface aerodynamic theory in Chapter 4, and 
the sensitivity of flutter speed to shape and modal parameters is computed using 
automatic differentiation (AD IFOR). Finite element modeling of wings is resorted 
to in Chapter 5 for structural modeling to generate the free vibration mode shapes 
required for aeroelastic analysis. The derivatives of flutter speed with respect to 
shape parameters are computed using a combination of central difference scheme 
and ADIFOR and the sensitivity to modal parameters are calculated using ADI- 
FOR. An equivalent plate model incorporating first-order shear deformation theory 
is used in Chapter 6 for structural modeling and the natural frequencies and flut- 
ter speeds obtained using the first-order shear deformation theory are compared 
with those obtained using the classical laminated plate theory. The sensitivity of 
natural frequencies to changes in shape parameters is obtained using ADIFOR. 
In order to demonstrate the use of shape sensitivities, an optimization problem of 
minimum weight design of a wing subject to flutter constraints is also presented. 
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CHAPTER 2 


SENSITIVITY OF AEROELASTIC RESPONSE 
OF A TYPICAL SECTION 

2.1 OVERVIEW 

In recent years, considerable effort has been made to integrate the aerody- 
namic, structural and control aspects of the design of an aircraft . Since the control 
and the structural dynamic behaviors can easily be expressed in the state-space 
form (i.e,, in terms of a set of first order ordinary differential equations in time), it 
is desirable that the unsteady aerodynamic airloads be also expressed in the same 
form. 

The state-space approach has the advantage that any system of differential 
equations can be represented by a set of first order ordinary differential equations 
of the form 

x = Ax + Bu 

with the output equations given by 

y = Cx + Du (2.1) 

where x are the aerodynamic state variables, u are the system inputs and y are 
the system outputs. If the unsteady aerodynamic behavior can be represented 
by state equations, then they can be easily coupled to the structural equations of 
motion and the resulting system can be examined for aeroelastic stability. 
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The unsteady aerodynamic loads on a typical section in response to an arbi- 
trary forcing can be represented through the use of indicial (step) response func- 
tions. By definition, an indicial function is the response to a disturbance that is 
applied instantaneously at time zero and held constant thereafter, i.e., a distur- 
bance given by a step function. A main advantage of the approach is that when 
the indicial response to a particular forcing mode is known, (e.g., that due to angle 
of attack or pitch rate) the cumulative response to an arbitrary forcing function 
can be obtained in the time-domain by means of Duhamel superposition. For in- 
compressible flow, the indicial lift response was first derived by Wagner [20] and 
is known exactly in terms of Bessel functions. However, since it is impractical 
in several applications to repetitively evaluate the Bessel functions, the Wagner 
function is usually approximated using sum of exponential functions. For example, 
Jones [21] used a two-pole exponential approximation to the Wagner function [20] 
given by 

<f>(S) = 1 - 0.165erp(— 0.04555) - 0.335exp(-0.35) (2.2) 

where 5 = 2 Vt/c, V is the freestream velocity, t is the time and c is the chord. 
The state equations describing the unsteady aerodynamic response can then be 
obtained by the application of Laplace transforms to these indicial functions. The 
resulting state equations are 


-0.01375(2f ) 2 — 0.3455(^f ) J \ r 2 


!!) + {? K/4(0 


with the output equation given by 


21/ , 2k f ci 1 

C N {t) = 2tt[0.006825( — f 0.10805( — )K 1 }+0.5a 3 /4(*) 

C C ( J ' 
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where Cm is the normal force coefficient and a 3 / 4 (f) is the quasisteady angle of 
attack at 3/4 chord. 


2.2 AERODYNAMIC MODEL 


In this analysis, the state-space representation given by Leishman and 
Nguyen [23] has been used to represent the compressible unsteady aerodynamics. 
They have represented the aerodynamic indicial response functions for compress- 
ible flow by up to three-pole approximations, the response consisting of two parts, 
one due to non-circulatory loading and the other due to circulatory loading. The 
indicial normal force and quarter chord pitching moment responses for a typical 
section to a step change in angle of attack a and a step change in pitch rate q can 
be written as [23] 


= ^Hs,M) + C N SM)4> c a {S,M) 

= -4^m(5,M) + C n .(M)^(S,M)( 0.25 - z«(M)) 
am 




-•i— 


(2.5) 


Cm(S) 7 i , . Cn 0 (M) c 

— T— = - (p qM [S,M) 


q J 16 

where <f>%, <^, <p£y<^, $ \m are exponential functions of S and M. Here, 


M is the Mach number, q = ac/V is the pitch rate, Cm is the normal force 
coefficient, Cm is the pitching moment coefficient about the quarter chord and 
C M a is the normal force curve slope. The superscripts C and I refer to circulatory 
and non-circulatory components of the indicial response functions. 'Note that ^ 
(where = \/l — M 2 is the compressibility factor) in [23] has been replaced by 
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CV a (M), so that experimental values of CjV a obtained as functions of Mach number 
can be used. 


The aerodynamic state equations have been shown by Leishman and Nguyen 
[23] to be given by 

x = Ax + B < 


a 

<7 


where 


A = diag[a\i 0,22 033 a 44 <255 066 077 ass ] 

,T 

11 1 i 111 1 11 n j 

B = 


1 110 110 0 

0.5 0.5 0 1 0 0 1 1 


The output equations are given by 


( 2 . 6 ) 


/ C * 

1 Cm 


Cx + D { a 1 

1 q J 


(2.7) 


where 


C = 
D = 


C11 Cl 2 Cl 3 Cl 4 
C21 C22 0 0 

4/M 1/M 

-1/M —7/12 M_ 


0 

C25 


0 

C26 


0 0 

C27 C28 


The nonzero terms of the a;/s and Ci/s are given in the Appendix. 


2.3 AEROELASTIC MODEL 


The aerodynamic equations in state-space form can be coupled to the struc- 
tural equations of motion of an airfoil section with bending and torsional degrees 
of freedom. The equations of motion for the airfoil section shown in Fig. 2.1 can 
be written as 

mh + Sgd + gu h + mu\h = Qh 


Sgh + Ig 9 + ggO + Igujgd — Qg 


( 2 . 8 ) 
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where m = fp,p(c/2) 2 is the mass per unit length, p is the mass ratio, p is the air 
density, Ig = m(c/2) 2 rj is the polar moment of inertia about the quarter chord 
per unit length, rg is the radius of gyration about elastic axis, So = m(c/ 2)xg 
is the static mass moment, xg is the nondimensional distance in semichords from 
elastic axis to center of mass, c is the chord length, h is the plunge displacement 
(positive downward), 0 is the pitch angle, u>h and ug are the bending and torsion 
frequencies respectively, gh and gg are the structural damping coefficients in plung- 
ing and pitching respectively and Qh and Qg are generalized aerodynamic forces 
in plunging and pitching respectively. 

By defining the states 


z\ 


h, z 2 = 0, 23 


h, 24 = 0 , 


( 2 . 9 ) 


the above equations can be written as 


where 


I 

0 



0 

-k 



M = 


m Sg 


9h 

0 ' 

_ Sg Ig 

g ^ 

_ 0 

90, 


Q — 


vribj\ 0 

0 IgUJg 


-L 

M 



-c 0 
0 c 2 



( 2 . 10 ) 


In order to couple the structural and aerodynamic equations, the input vector 


can be expressed in terms of the 2 states as given below 



'0 1 l/V 

0 

l*/" 

0 0 0 

c/V_ 




*% 

23 

24 


( 2 . 11 ) 
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2.4 SENSITIVITY EQUATIONS 

The aeroelastic equations obtained as a set of first order ODEs is of the form 

IP] w = [Q] W (2.14) 

which could be written as 

w = [E] w (2.15) 

where [ E ] = [P] -1 [Q] 

Since the flutter speed is determined by solving the eigenvalue problem, the 
derivatives of the eigenvalues have to be calculated for obtaining the flutter speed 
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derivatives. A number of papers [76-80] have been published on calculation of the 
derivatives of eigenvalues and eigenvectors of real, symmetric matrices as well as 
general matrices. 

The derivative of the ith eigenvalue for the eigenproblem in equation (2.15) 
with respect to the flutter speed is given by [81] 


ax‘ { 4 } 

dV f {ejFfej.} 


(2.16) 


where {e|} and (e)} axe the kh left and right eigenvectors respectively, and -QP is 
calculated analytically by differentiating the elements of the matrix with respect 


to Vf. 


Similarly, the derivative of the ith. eigenvalue with respect to any parameter 
p is given by 

(2.17) 


dp 


dp {ejF{4} 


can be conveniently written as 




(2.18) 


and can be computed analytically, where 




12.19) 


The ith eigenvalue \ l is a function of the speed V and the parameter p which 
are varied. That is, 

A* = V(V,p) (2.20) 
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Therefore, 



Thus we have, 

0 = Real( — - )dV f + Real( — )dp (2.23) 

oVf op 

The analytical derivative of the flutter speed with respect to parameter p is 
then given by 

= _ Mfl (2.24) 

dp *<■'(!#) 

The [E] matrix is composed of mass, stiffness and aerodynamic matrices. 
Obtaining the analytical derivatives of the mass, stiffness and aerodynamic terms 
with respect to any parameter p is straightforward for the typical section. This 
is done by computing analytically the derivatives of those terms that are explicit 
functions of the desired parameters. 


2.5 AEROELASTIC ANALYSIS 

An aeroelastic formulation in the time-domain leads to a set of first-order 
ordinary differential equations (ODEs) which can be solved by time-integration or 
as an eigenvalue problem to examine the stability of the system. 
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2.5.1 EIGENVALUE SOLUTION OF THE AEROELASTIC EQUATIONS 


The flutter characteristics of the airfoil are found by calculating the complex 
eigenvalues A f. = erf. + ia>k at various values of free stream velocity. Flutter occurs 
at the lowest speed for which any oq- becomes positive. 

The flutter speed was determined for the following case, the results for which 
have been presented by Leishman and Crouse [82]. The parameters used are /j, = 


100, xq = 0.25, rg = 0.5, u>h = 10 rad/s., ujg = 50 rad/s., ah — —0.5,6 = 


5 in., CN a = 14.65, x ac = 0.286, M — 0.85. Flutter was found to occur at 92.34 
ft/s., i.e., a non-dimensional speed of V /bu>g=^AZ which agrees well with the value 
of V /bug—- 4 .4 reported in [82]. 

The damping ratio for each of the aeroelastic modes is given by 


v ' 7 k + 

A plot of the variation of the damping ratio ( with non-dimensional speed V/bujg 
is given in Fig. 2.2. Fig. 2.3 shows the variation of flutter speed predicted by this 
method for different values of u>h/ug. 



2.5.2 TIME INTEGRATION OF THE AEROELASTIC EQUATIONS 


Several numerical integration techniques [83-85] exist for step-by-step solu- 
tion of these equations. The aeroelastic equations are integrated with respect to 
time using a time-integration scheme. The Wilson-0 method [85] was used for this 
purpose. A set of first order ODEs can be represented as 


MM = [S\{*} 


(2.26) 
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In the Wilson-0 method, it is assumed that the variation of velocity from time t 
to t + 9 At, where 0 > 1.37, is linear. At time (t + 0A t), then 

Q A J. 

{ut+e&t} = {«*} + -^-Wt+e&t + u' t } (2.27) 

Then equation (2.26) becomes 

[(«] - [5]^] M +tol } = tsu«. + ^<4J (2.28) 

Using the starting values of {u<} and {u' f } at time t, {u' t+dAt } is computed from 
equation (2.28). The vector {u t +BAt} is then calculated from equation (2.27). The 
new values of {ut+s&t} and {u' t+ $ At } are then used in equation (2.28) to update 
the {u' t+9At } vector. The step-by-step integration of the equations is done in this 
manner with respect to time by repeating the above process. The amplitudes of 
displacements are then monitored as time progresses. 

Flutter for the typical section (see section 2.5.1) was found to occur at 92.34 
ft/s., i.e., a non-dimensional speed of V/bojg=4:A3. A time-integration of the 
first order ODEs representing the aeroelastic system was done using the Wilson-0 
method. The plunge and pitch amplitudes of motion are plotted with respect to 
time in Figs. 2.4 and 2.5 respectively, at three different non-dimensional speeds, 
including the flutter speed. It can be seen that at speeds below the flutter speed, 
the oscillations that are set in due to any initial disturbance given to the airfoil die 
out as time progresses, whereas, at speeds above the flutter speed, the displace- 
ment amplitudes increase with time, leading to instability. At the flutter speed, 
the oscillations are able to maintain a constant amplitude, denoting a neutrally 
stable condition. 


20 



2.6 SENSITIVITY RESULTS 


The sensitivity of flutter speed with various parameters namely yu, xg, rg, u>/,. 
and ug was calculated by both analytical and finite difference methods. In the an- 
alytical method, the derivatives of the [E] matrix (see section 2.4) with respect 
to the above mentioned parameters were calculated analytically. The finite dif- 
ference derivatives were calculated for step sizes of 1%, 0.1% and 0.01%. The 
parameters were perturbed one at. a time using these step sizes and the flutter 
speed recomputed. A forward difference scheme was then applied to compute the 
derivatives. It can be seen from the results shown in Table 2.1 that the forward 
difference derivatives obtained using a step size of 0.01% have good agreement with 
the analytical values. 

Figs. 2.6-2.10 show the variation of flutter speed obtained by eigenanal- 
ysis with respect to various parameters. In each case, the sensitivity derivative 
computed at the baseline configuration is also shown. The sensitivity derivative 
computed forms a tangent to the flutter speed curve at the baseline, and in some 
cases, the linear approximation given by the sensitivity analysis approximates the 
flutter speed curve very closely, over the entire range of the parameter that is 
varied. 


21 




Nondimensional speed, V/b(D e 

Fig. 2.2 Variation of damping ratio, £ with 
nondimensional speed, V/bco a 








Plunge displacement 




Pitch displacement 





Flutter speed 
Sensitivity analysis 



r speed vs Mass ratio (M=0.85) 






Flutte 



Static unbalance, x Q 

0 


Fig. 2.7 Flutter speed vs Static unbalance (M=0.85) 
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Fig. 2.9 Flutter speed vs Bending frequency (M=0.85) 
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Table 2.1 Sensitivity of flutter speed with respect to various 
parameters at M=0.85 (Logarithmic derivatives) 

(jj. = 100, xg = 0.25, rg = 0.5, Uh = 10 rad/sec., ujg = 50 rad/sec., 
6 = 5 in., ah = —0.5, V/ = 92.34 ft/sec) 


Parameter 

Analytic 

derivative 

Finite Difference 
derivative 

1%“ 

0.1% 

0.01% 

M 

0.42556 

0.42016 

0.42499 

0.42548 

XB 

-0.68661 

-0.67349 

-0.68525 

-0.68644 

re 

1.29638 

1.28080 

1.29474 

1.29615 

u>k 

-0.47110 

-0.46741 

-0.47071 

-0.47104 

we 

1.47113 

1.45548 

1.46948 

1.47089 


a indicates step size 
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CHAPTER 3 


SENSITIVITY OF AEROELASTIC RESPONSE OF A WING 
USING TIME DOMAIN APPROACH 

3.1 OVERVIEW 

In the previous chapter, a state-space model for the unsteady aerodynamics 
was used for aeroelastic analysis of a typical section and the sensitivity of flutter 
response to airfoil parameters were computed analytically. An aircraft during flight 
is subjected to unsteady airloads and the unsteady lift forces and moments at any 
spanwise station of the wing can be represented by the state-space formulation for 
the typical section presented earlier. A strip-theory formulation is thus developed 
in this chapter, to represent the unsteady aerodynamic forces on the wing. The 
wing structure is modeled using classical plate theory and is based on a Rayleigh- 
Ritz formulation using Chebyshev polynomials for the wing displacements. The 
structural equations are coupled with the aerodynamic state equations and the 
resulting aeroelastic equations are solved as an eigenvalue problem to examine the 
onset of aeroelastic instability. The equations are also integrated with respect to 
time to observe the aeroelastic phenomena in real time, at speeds close to the 
critical speed. The structural and aerodynamic characteristics of the wing are 
functions of its shape and hence the flutter response is sensitive to changes in 
shape parameters. In this chapter, the sensitivity of aeroelastic response of the 
wing with respect to the shape parameters, namely, aspect ratio, area, taper ratio 
and sweep angle are calculated analytically. 
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3.2 STRUCTURAL MODEL FOR THE WING 


The structural formulation is based on a Ritz solution technique using the 
energy functionals for a laminated plate which includes the bending and stretch- 
ing of the reference surface. The planform geometry can be represented by any 
generally tapered skewed configuration. The original rectangular (x,y) coordinate 
system and the transformed (77, £) coordinate system of the wing are shown in 
Fig. 3.1. The x-y plane is the mid-plane of the wing and the z axis is normal 
to the wing. For an unswept wing the fiber angle is measured counterclockwise 
from the positive y axis. As the wing is swept, the fiber angle is also rotated 
correspondingly. 

In the Rayleigh- Ritz formulation, Chebyshev polynomials T; are used to rep- 
resent the displacements at any point on the wing [86]. The Chebyshev polynomials 
are given by 

W) = 1 

TiW = V> (3.1) 

Ti(<P) = 2m-i -Tt-i - 1 < < 1 

The displacements are expressed in terms of the Chebyshev polynomials as shown 

I J 

v(v,o = 'E.'E R n T iWm) 

i= 0 /== 0 
K L 

V( n , 0 = EE SMvWO (3.2) 

t=o /=o 
M N 

w(, 7,0 = EE P mn T m ( V )T„(B -1 <>),«<! 

m=0 7i=0 
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It has been shown by Singhvi and Kapania [86] that for free vibrations of 
the laminated composite wing (i.e., in the absence of aerodynamic forces) the 
equations of motion can be derived using classical plate theory in the form 

[«]{«} + [£]{«} = 0 (3.3) 

where [A'] and [M] are the stiffness and mass matrices, respectively. The vector 
{g} is defined as 

{<?} = ( AoOj Aoi j •••Rij] Sqq, Sqi , -PqOi Toi, ... Pmn ) T (3.4) 

Linear and rotational springs of large magnitude are placed at the wing root to 
approximately satisfy the clamped boundary conditions. The stiffness matrix for 
the plate alone (i.e., excluding the springs) is 

[A-] = //VlW g £ [TPH-Wf (3.5) 

where [ B ] is the matrix whose elements consist of the partial derivatives of the 
Chebyshev polynomials with respect to the natural coordinates r) and £ and is 
defined by 

(*',) = [B]{«} (3-6) 

where 

J = (u, v n V£ W n W V r, Wtf w v $) 

The [T] in equation(3.5) is the transformation matrix that relates the strain and 
curvature vector in the ( x,y ) coordinate system to the strain and curvature vector 
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in the (?y,£) coordinate system and ./ is the Jacobian of the transformation. The 
strain transformation is given by 




(3.7) 


where 


('U;£ Vy Uy + V x tOxx ^ yy X y ) 


The details of the [T] and [. B ] matrices and J are given in [87]. A typical element 
of the mass matrix [M] is given by 


Mij = Pm t J f T k ( v )T,(i)T 0 (r,)T r mJ\d v d(, (3.8) 

The coefficients Rij and Ski in {q} corresponding to the inplane displace- 
ments in equation(3.3) are condensed out using static condensation to the form 


[M){P mn } + [K]{P mn } = 0 


(3.9) 


where [M] is the mass matrix and [K] is the stiffness matrix of order (M-f-l)x(iV-j-l) 
(see equation(3.2)) with generalized coefficients {P mn }- In the present work, a 
value of 5 is chosen for both M and N. 


3.3 AEROELASTIC MODEL FOR THE WING 

The aerodynamic state space model which was used for the aeroelastic anal- 
ysis of a typical section is extended to represent the unsteady aerodynamic forces 
acting on a wing. The lift and moment forces on a typical section acting at the 
quarter chord are given by equation(2.10) as 

L = | pV 2 cC N 
M = \ P V 2 c 2 C m 
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When extending this compressible aerodynamic theory to a finite span wing, 
the lift forces are assumed to be distributed along the quarter chord line (reference 
line) and the moments act about the reference line. Since the lift and moment 
forces are non-conservative forces, using the principle of virtual work, we get 

SW nc = f —LShdy + j MSddy (3.10) 

Jo Jo 

where l is the length of the quarter chord line, 6h and 86 are virtual displacements 
and y is the coordinate along the reference line. 

The displacement at any location y is given by 

h(y) = w (??,£) (3.11) 

where rj and £ are the natural coordinates corresponding to the (x, y) coordinates 
of the point at distance y from the origin. 

The rotation about the reference line (positive wing leading edge up) is given 
by 

6(y) — w rx cosA — w,g sini\ (3.12) 

To facilitate numerical integration using Gaussian quadrature, the limits of in- 
tegration along the reference line are transformed in the range of —1 to 1 by 
y == 1(1 + ■0)/2, where —1 < ip ^ 1 


Substituting the expressions for the lift and moment on the wing and the 
wing deflection, we have 
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and 


f M89dy 

Jo 

1 f 1 

= - / MSOdJ; 

2 J-i 

/I m n r «i 

= 5</ ! >EE / c2 + Pvi 

i =0 ?=0 •* ^ 


% 0 


0 H 


ij J 


n 


) w 2 o# 


8Pij 

(3.13) 


where m and n represent the order of the Chebyshev polynomial used in the 
displacement function. The row vectors [Ci p ] and [C 2 P ] are the elements of the [C] 
matrix (see section 2.2) where p — 1,2, . . . ,8. The row vectors [D\ p '\ and [T* 2 p'] 
are the elements of the matrix given by 


[D] = 


‘ 4/M 

1/M 

[0 1 1/7 0 1 

_ -1/M 

-7/12 M. 

o 

o 

o 


(3.14) 


(3.15) 


where p' = 1, ... ,4. 

The variables and W 2 tj in equation(3.13) are given by 

w Ul = TiWTjiO 

u> 2 ,.y = cosA(Ti )V Tjrj, x ) sinA{Ti tr JTjVj^y ) 

where A is the sweep angle. 

[Bn] in equation(3.13) is a matrix of order 2xN where N = (m + l)(n + 1). 
A typical column of the matrix is given by 

(3.16) 


ff = I TiTj 

\ cosA(Ti <r ,Tjr], x +TiTj^i, x ) — sinA(Ti, n Tjr), y +T t Tj^, y ) 


where k — 1,2, . .. . , N . 

The column vector {r} in equation(3.13) is the vector of aerodynamic state 
variables and {P t] Pij} T is the vector of generalized displacements. 
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Using equations (3.11), (3.12) and (3.13), equation(3.10) can be written as 

N 

SW nc = QitPi 1 N = (m + l)(n + l) (3.17) 

i=i 

where 

Qi = [C']{x} + (d; Dilj ^ ) (3.18) 

The aerodynamic state equations(2.6) for a typical section perpendicular to 
the quarter chord line were in the form 


{i} = [A]{z} + [£]{u} 


An integration of these state equations along the quarter chord line to consider 
the effect of finite span yields 






H. 


ij o 

0 H \ 


*j J 




*3 


d'tf) 


= [A']{x} + [B[ B'o}< 


u 




(3.19) 


Equations(3.9) can be written as a set of first order ODE’s in {Pij} and 
{Pij} which will be represented by {p l } and {<?,}, respectively. It can be coupled 
with equations(3.18) and (3.19) to generate the aeroelastic equations of the wing 
in the form 


■/ 

0 

O' 

{ & } 

0 

1 

0 ' 

0 

M 

0 

< > = 

D\ — K 

a 

0 

0 

I 

\ x J 

B[ b[ 

A! 



(3.20) 

Since a Chebyshev polynomial of order 5 is chosen for the displacement function 
in q and £, we have 36 generalized coefficients {p;}, their 36 time derivatives {qi} 
and the 8 aerodynamic state variables {a;}. The stability of this system can be 
determined by solving this 80x80 eigenvalue problem. 
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3.4 NUMERICAL EXAMPLES 


The natural frequencies and flutter speeds predicted by this aeroelastic state- 
space formulation are compared with previously published results. The aeroelastic 
equations are solved as an eigenvalue problem and also using a time-integration 
scheme. 

3.4.1 EIGENVALUE SOLUTION OF THE AEROELASTIC EQUATIONS 

Before performing the sensitivity calculations of the flutter speed of the wing 
with respect to shape parameters, comparison of natural frequencies and predicted 
flutter speeds was made with other results. The first three natural frequencies of 
an unswept wing and its flutter speed in subsonic flow were compared with results 
reported by Landsberger and Dugundji [88] for different laminate sequences in 
Table 3.1 and Table 3.2, The material used for the wing is Hercules Graphite 
Epoxy (AS1/3501-6) with properties: E\ — 98 x 10 9 Pa, Eo. = 7.9 x 10 9 Pa, 
V \2 = 0.28, G \2 — 5-6 x 10 9 Pa and p — 1520 kg/m 3 . The thickness of each ply 
is 0.134 x 10 ” 3 m. The flutter data used for comparison [88] are the experimental 
results from the wind tunnel tests performed in the MIT Acoustic wind tunnel. 
The results agree fairly well. 

3.4.2 TIME INTEGRATION OF THE AEROELASTIC EQUATIONS 

The aeroelastic equations are integrated with respect to time to examine the 
effect of time-dependent airloads on the vibration characteristics of the structure. 
This gives insight into the type of response that can be expected from a flexible 
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structure operating at conditions on the verge of instability. One of the features 
of time-integration is that one can observe the aeroelastic response in real time 
(important for studying fatigue) at different airspeeds (he., at different freestream 
velocities), than just determining the critical speed at the onset of instability. 

The aeroelastic response of the wing shown in Fig. 3.2 is examined. The 
wing skins are made of 0 deg. laminated Graphite/Epoxy (T300/iV5208) with 
the following material properties: E\ = 181 x 10 9 Pa, Ei = 10.3 x 10 9 Pa, 
1/12 = 0.28, (?i 2 = 7.17 x IQ 9 Pa and p = 1600 kg/m 3 . The aeroelastic analysis 
is carried out at Mach number of 0.7 and pair = 0.6 kg/m 3 . In order to observe 
the aeroelastic phenomena in real time, one of the coefficients of the displacement 
function was perturbed and the system of equations was integrated with respect 
to time using the Wilson-0 method, described in Section 2.5.2. The wing tip 
displacement is plotted as a function of time at different speeds in Figs. 3. 3-3. 9. 
The tip displacement of the 5° swept wing at the divergence speed is shown in 
Fig. 3.3. It can be seen that the displacement approaches a constant amplitude 
at this speed and does not decrease to the zero mean value. Above this speed, the 
displacement increases with time as shown in Fig. 3.4. This is characteristic of 
divergence which is a static instability, and hence of a non-oscillatory nature. The 
small wiggles in the response are due to the dynamic effects. For the 20° swept 
wing, the displacement slowly approaches a constant value at the divergence speed 
as shown in Fig. 3.5. The oscillatory nature is due to the fact that the flutter 
and divergence speeds are close enough. It can be seen that the oscillations die 
down since the dynamic instability has not yet set in. Fig. 3.6 shows the flutter 
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condition for the 20° swept wing. Since the wing has already diverged at a lower 
speed, we see the constant amplitude oscillations about a diverging mean position. 
This is the condition when a static mode has become unstable and the oscillatory 
mode is neutrally stable. Figs. 3.7, 3.8 and 3.9 show the tip displacement of the 
30° swept wing below the flutter speed, at the flutter speed and above the flutter 
speed, respectively. The oscillations are damped out at low speeds, but at the 
critical speed the structure is able to maintain self-sustained oscillations, and at 
higher speeds, the amplitudes of oscillations rapidly increase with time. 

3.5 SENSITIVITY RESULTS 

The aeroelastic equations (3.20) are a set of first-order ordinary differential 
equations and the procedure for analytically determining the sensitivity of the flut- 
ter speed is described in section 2.4. The expressions for the analytical derivatives 
of the [A'] and [M] matrices (equations(3.5) and (3.8)) are given in [87]. Since the 
reduced stiffness matrix [A'] is obtained from [A'] by static condensation the ana- 
lytical derivative is obtained by a succession of differentiations using the chain 
rule. The derivatives of the aerodynamic terms are obtained by taking the analyt- 
ical derivatives of those terms that are explicit functions of the shape parameters, 
given in Appendix A. 

Sensitivity analysis of the flutter speed of the wing with respect to shape 
parameters is carried out for the wing shown in Fig. 3.2. The wing skins are 
made of 0° laminated Graphite/Epoxy (T300/iV5208) with the following material 
properties: E\ = 181 x IQ 9 Pa, E% = 10.3 x 10 9 Pa, v \2 = 0.28, G 12 — 7.17 x 
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10 9 Pa and p = 1600 kg/m 3 . The flutter speed calculated using this strip-theory 
formulation is compared with the flutter speed obtained using a lifting-surface 
unsteady aerodynamic theory (FAST) in Fig. 3.10 for the wing at Mach 0.7 and a 
range of values of aspect ratio. The results agree fairly well, with the strip-theory 
predicting higher values for the flutter speed. The critical airspeed of the wing is 
shown in Fig. 3.11 as a function of the quarter-chord sweep angle. As seen from 
the graph, divergence (zero frequency flutter) instability is critical upto a sweep 
angle of about 20° and for higher sweep angles, the flutter mode is the unstable 
mode. The shape sensitivity derivatives of the divergence speeds and flutter speeds 
of the wing at Mach 0.7 and p a ,> = 0.6 kg/m 3 are computed analytically. 

The critical speeds of the wing obtained by perturbing one shape parameter 
at a time from the baseline configuration are shown in Figs. 3.12-3.19. The 
prediction of critical speed by analytical sensitivity calculations is also superposed. 
The sensitivity derivative obtained forms a tangent to the critical speed curve at 
the value of the shape parameter at which it is computed. 

It is observed from Fig. 3.12 and Fig. 3.16 that the divergence and flutter 
speeds, respectively, drop as the aspect ratio increases. Physically, as the wing 
is made more slender, the instability sets in at a lower speed. As the area of 
the wing is increased, the divergence and flutter speeds are found to decrease as 
shown in Fig. 3.13 and Fig. 3.17, respectively. The greater the area of the wing, 
the greater the unsteady aerodynamic forces acting on it, causing the aeroelastic 
instability to occur at a lower speed. The analytical sensitivity gives a good linear 
approximation of the critical speeds over the range the aspect ratio and area of 
the wing are varied. 
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The divergence speed of the wing is not affected much as the taper ratio is 
increased from 0.4 to 0.6 as seen from Fig. 3.14. The flutter speed of the wing 
decreases as the taper ratio is increased as shown in Fig. 3.18. The flutter speed 
curve can be approximated linearly as indicated by the sensitivity prediction. 

It can be seen from Fig. 3.15 that the divergence speed increases as the wing 
is swept backwards. The flutter speed does not vary much with increasing sweep 
angle about the 15° swept configuration as seen from Fig. 3.19. By performing 
one sensitivity calculation at the baseline analytically, this method gives a linear 
approximation to the critical speeds of the wing for changes in the wing shape 
parameters about the baseline. This information is useful for preliminary design 
purposes, as it avoids the necessity of a reanalysis for small changes in any of the 
shape parameters. 
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Fig. 3.5 Displacement vs Time at V=92.64 m/s 

(AR=10, Area=20m 2 , TR=0.5, Sweep=2(f) 
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Fig. 3.6 Displacement vs Time at V=95.90 

(AR=10, Area=20m 2 , TR=0.5, Sweep=2(f ) 
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Fig. 3.7 Displacement vs Time at V=90.0 m/s 

(AR=10, Area=20m 2 , TR=0.5, Sweep=3(f ) 
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Fig. 3.8 Displacement vs Time at V=98.87 m/s 

(. AR=10 , Area=20m 2 , TR=0.5, Sweep=3(f) 
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Fig. 3.9 Displacement vs Time at V=100.0 m/s 

(AR=10, Area=20m 2 , TR=0.5, Sweep=3(f) 
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Fig. 3.10 Flutter speed from strip-theory (state-space aerodynamics) 
compared to lifting-surface theory (FAST) 

(M=:0.7, Area=20m 2 , TR=0.5, Sweep=Cf) 
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Fig. 3.11 Critical speed vs Sweep Angle 

(AR=10, Area=20 rrf, TR=0.5) 
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Fig. 3.13 Divergence speed vs Area (M=0.7) 

(AR=10, Area=20m 2 , TR=0.5, Sweep=0°) 
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Fig. 3.18 Flutter speed vs Taper ratio (M=0.7) 

(AR=10, Area=20m 2 , TR=0.5, Sweep=15°) 
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Fig. 3.19 Flutter speed vs Sweep angle (M=0.7) 

(AR=10, Area=20m 2 , TR=0.5, Sweep=15°) 





Table 3.1 Comparison of natural frequencies of an unswept wing 
(Area = 0.02318 m 2 , Aspect ratio — 4.0132, Taper ratio — 1.0) 
for different laminate sequences 


Laminate 

sequence 

Natural frequencies (Hz.) 

Present 

Landsberger and Dugundji [88] 

First 

Second 

Third 

First 

Second 

Third 

[0 2 /90], 

11.03 

39.30 

69.06 

10.8 

39 

67 

[152/0], 

8.86 

42.62 

63.25 

8.5 

48 

58 

[±15/0], 

10.12 

48.9 

.64.94 

9.9 

50 

63 

[+30 2 /0], 

6.21 

37.57 

57.78 

6.0 

41 

60 

[±30/0], 

7.73 

48.76 

64.42 

7.8 

50 

65 


Table 3.2 Comparison of flutter speed of an unswept wing 
(Area = 0.02318 m 2 , Aspect ratio == 4.0132, Taper ratio = 1.0) 

in subsonic flow 


Laminate 

sequence 

Flutter speed (m/s) 

Present 

Experimental [88] 

[0 2 /90], 

24.9 

26 

[15 2 /0]s 

23.2 

25 

[±15/0], 

28.1 

28 

[T15/0], 

18.9 

21 

[+30 2 /0], 

28.8 

29 
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CHAPTER 4 


SENSITIVITY OF AEROELASTIC RESPONSE OF A WING 
USING FREQUENCY DOMAIN APPROACH 

4.1 OVERVIEW 

The problem of flutter instability was studied in Chapter 3 using a state-space 
unsteady aerodynamic representation and a Rayleigh-Ritz structural formulation. 
The sensitivity of flutter speed to shape parameters of the wing were obtained 
analytically. In this Chapter, the same Rayleigh-Ritz structural formulation is 
adopted, but the aeroelastic analysis is carried out using lifting-surface unsteady 
aerodynamics using modules from a system of programs [89] called FAST (Flutter 
Analysis System). Flutter speeds are obtained using a V-g type of solution. The 
derivatives of the flutter speed with respect to shape parameters of the wing are 
calculated using automatic differentiation (ADIFOR). In order to determine the 
importance of a particular mode to flutter, derivatives of the flutter speed with re- 
spect to natural frequencies, generalized mass and generalized aerodynamic forces 
are computed using ADIFOR. 

Automatic differentiation is emerging as a valuable tool for sensitivity cal- 
culations. ADIFOR, GRESS, PADRE-2, Power Calculus and ODYSEE are some 
of the automatic differentiation packages [90] developed for differentiating FOR- 
TRAN77 codes. ADIFOR (Automatic Differentiation of FORtran) is a joint ven- 
ture of Rice University and Argonne National Laboratory. ADIFOR processes a 
given Fortran code and generates a Fortran code for computing the derivatives 
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of the desired output variables with respect to the independent variables by ap- 
plying the chain rule of differentiation. Differentiating large codes by hand is 
cumbersome, divided differences are dependent on choice of step size and sym- 
bolic programs may be infeasible for large codes. Automatic differentiation, on 
the other hand, can handle codes of arbitrary size and produce exact derivatives 
of the discrete system with no truncation error. 

4.2 AEROELASTIC MODEL 

The equivalent plate model described in section (3.2) based on a Rayleigh- 
Ritz formulation using Chebyshev polynomials for the wing displacements is used 
for structural modeling. Flutter calculations are performed using FAST [89], a 
system of programs based on the subsonic kernel function lifting-surface aerody- 
namic theory. A free vibration analysis of the wing is performed and the vibration 
modes from this analysis are fed into the modes processing module from FAST. 
The subsonic kernel function matrix program then solves the subsonic downwash 
integral equation for the oscillating planar wing lifting surface. The generalized 
force module from FAST then computes the aerodynamic forces from the subsonic 
kernel matrices. The flutter speed of the wing is obtained using a V-g type of 
solution. The generalized aerodynamic forces are determined for a specific Mach 
number and for a range of values of the reduced frequency for the specified down- 
wash distribution. These values are then interpolated to get aerodynamic forces 
for closely spaced values of the reduced frequency. The flutter equation solved by 
the program is 

n 

[w 2 -uf( 1 + ig))Miqi + ^ A^qj — 0, i = 1, n (4.1) 

i=i 
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where w is the vibration frequency, «; is the frequency of the zth natural vibration 
mode, Mi is the generalized mass associated with the zth natural vibration mode. 
g is the incremental damping, Aij are the generalized aerodynamic forces resulting 
from the pressure induced by the yth mode acting through the displacements of 
the ith mode and qi is the ith component of the flutter eigenvector. 

In terms of the non-dimensional generalized aerodynamic forces Aij , the 
above equation can be written as an eigenvalue problem in the form, 

n 

Cjj - SijQ,)qj == 0, i = 1, n (4.2) 

7=1 

where the eigenvalue Q is given by 

2 

Q= (~V Z ) + ( 4 - 3 ) 


and 



( V T P b l A ij 

\u n) [ 2 Mi 


+ Sijk 2 


(4.4) 


where ui a is a reference frequency, b 0 is the reference length, usually root semichord, 
k is the reduced frequency, p is the air density and V is the airspeed. 

The above eigenvalue problem is solved for a range of values of the reduced 
frequency and monitored for crossings on the V-g plane when the incremental 
damping g goes to zero. At each of these crossings, the values of the airspeed and 
the frequency are noted. The critical flutter speed is the lowest speed at which the 
damping of the structure goes to zero, 
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4.3 SENSITIVITY ANALYSIS USING AUTOMATIC DIFFERENTI- 
ATION 

The sensitivity of flutter speed of the wing to shape parameters namely 
aspect ratio, area, taper ratio and sweep angle is computed using ADIFOR. In 
order to examine the importance of a particular mode to flutter, the sensitivity 
of flutter speed with respect to the modal parameters, namely, natural frequency, 
generalized mass and generalized aerodynamic forces is calculated using ADIFOR. 

4.3.1 AUTOMATIC DIFFERENTIATION 

Once the input variables and output variables of interest are identified, AD- 
IFOR can be used to differentiate FORTRAN codes. ADIFOR propagates deriva- 
tive values rather than formulas as in an analytic derivation. Only those inter- 
mediate variables that are functions of the input variables and are required in the 
computation of the output variables are identified and differentiated by ADIFOR. 

A small example of an ADIFOR-generated FORTRAN code to compute the 
derivatives from a simple program main.f with one call made to the subroutine 
func.f is given in Appendix B. ADIFOR generates the code g_func.f which has 
statements to compute the derivatives of the output variable y with respect to the 
input variable array r(3). 

4.3.2 SENSITIVITY TO SHAPE PARAMETERS 

For all the results presented, the wing shown in Fig. 3.2 is used with the 
wing skins made of 0 deg. laminated Graphite/Epoxy (T300/iV’5208) with the 
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following material properties: E\ — 181 x 10 9 Pa, E 2 — 10.3 x 10 9 Pa, v \ 2 = 0.28, 
G 12 — 7.17 x 10 9 Pa and p — 1600 kg/m 3 . 

The sensitivity derivatives of the flutter speed with respect to shape parame- 
ters are calculated using ADIFOR. The sensitivity results and the results from re- 
analysis are shown in Figs. 4. 1-4.4 for the wing at M = 0.6 and p a i r — 0.8 kg/m 3 . 
The decrease in flutter speed due to increasing aspect ratio is seen to be gradual 
from AR = 8 to AR = 10 in Fig. 4.1. Then, a more sharper drop is flutter speed 
is seen as the AR is increased to 12. The sensitivity derivative however predicts 
the trend seen at the baseline configuration. The flutter speed decreases with in- 
crease in area and taper ratio as shown in Figs. 4.2 and 4.3. Increasing the sweep 
angle has the effect of increasing the flutter speed about the 30° swept-back con- 
figuration as shown in Fig. 4.4. By performing one sensitivity calculation at the 
baseline, this method gives a linear approximation to the flutter speeds of the wing 
for changes in the wing shape parameters about the baseline. This information is 
useful for preliminary design purposes, as it avoids the necessity of a reanalysis for 
small changes in any of the shape parameters. ^ 

4.3.3 SENSITIVITY TO MODAL PARAMETERS 

Complex wing structures are often modeled with a large number of degrees of 
freedom and a vibration analysis yields a large number of free vibration modes. A 
certain number .of these modes have to be fed into FAST to generate the generalized 
aerodynamic forces required for flutter analysis. Some of these inodes do not 
actively participate in the flutter and using several modes to determine flutter 
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instability leads to larger computational cost. In situations, where the natural 
frequencies and mode shapes of the wing are measured from experiments, one has 
large amount of modal data and has to determine the number of modes that are 
required for reasonably estimating the flutter speed. One could make a judicious 
choice of the number of modes that are required for flutter analysis, based on the 
sensitivity derivatives of the flutter speed with respect to the modal parameters of 
the wing. 

The sensitivity of flutter speed to modal parameters were obtained using 
ADIFOR. The derivatives of the flutter speed with respect to natural frequencies 
of the wing are shown in Fig. 4.5. Nine modes were used for the flutter analysis. 
It can be seen that if the third natural frequency is increased by 1 rad/s, then 
the flutter speed increases by 1,5 m/s. It is also seen that higher modes do not 
contribute much to the flutter. The derivatives of the flutter speed with respect 
to generalized mass are given in Fig. 4.6 and the flutter speed sensitivity to the 
real and imaginary parts of the generalized aerodynamic forces (GAF) is given in 
Fig. 4.7. The (i,j) term in Fig. 4.7 stands for the non-dimensional generalized 
aerodynamic force resulting from the pressure induced by the jth mode acting 
through the displacements of the ith mode. The sensitivities of flutter speed with 
respect to modal parameters give an estimate of the importance of a particular 
mode to flutter. 

In order to obtain a quantitative estimate of the contribution of a particular 
mode to flutter, one can use the sensitivity information which gives the variation 
of flutter speed with the natural frequency, generalized mass and real and imagi- 
nary parts of the generalized aerodynamic forces associated with that mode. It is 
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possible to construct a logarithmic derivative and sum up the absolute values of 
these derivatives, which reflects the change in flutter speed due to changes in modal 
parameters associated with a particular mode. Having obtained all the derivatives 
one could construct a parameter P; for each mode i given by 

p Vi dV f Mi dV f A T Re(Qij) dV f Im(Qy) dV f 

‘ V f dcoi + V f dMi^l V f dRe(Qij) + V f dIm(Q tJ ) 

(4.5) 

where Vf is the flutter speed, u>, is the natural frequency, Mi is the generalized 
mass and Re(Qij) and Im(Qij) are the real and imaginary parts of the generalized 
aerodynamic forces for an aeroelastic analysis using N vibration modes. The values 
of this parameter Pi for the nine modes used in the analysis are given in Table 
4.1. It can be seen that the higher modes do not contribute significantly to the 
flutter phenomenon. Based on this information, flutter analysis of the same wing 
was carried out using fewer modes (5 modes) and the flutter speed obtained was" 
100.338 m/s compared to a value of 100.404 m/s using 9 modes, which is not a 
significant drop in flutter speed. The differences in the values of the derivatives 
of the flutter speed computed using fewer modes are found to be insignificant 
compared to the derivative values using more number of modes. 
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Fig. 4.1 Flutter speed vs Aspect ratio (M=0.6) 

from FAST 

( AR=10 , Area=20m 2 , TR=0.5, Sweep=3CP) 
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Fig 4.2 Flutter speed Vs Area (M=0.6) 

from FAST 

(AR=10, Area=20m 2 , TR=0.5, Sweep=3(T) 
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Fig. 4.3 Flutter speed vs Taper ratio (M=0.6) 

from FAST 

(AR=10, Area=20m 2 , TR=0.5, Sweep=3CP) 
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Fig. 4.4 Flutter speed vs Sweep angle (M=0.6) 

from FAST 

(AR=10, Area=20m 2 , TR=0.5, Sweep-3CP) 
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Fig. 4.6 Flutter speed sensitivity to generalized mass (M=0.6) 

in (m/s)/(kg) 

(AR=10, Area=20m 2 , TR=0.5, Sweep=3(T) 
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Fig. 4.7 Flutter speed sensitivity to nondimensional generalized 
aerodynamic forces (GAF) (M=0.6) in (m/s) 

(AR=10, Area=20m 2 , TR=0.5, Sweep=3L f) 
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Table 4.1 Sum of the absolute values of the logarithmic derivatives 

(Parameter Pf for each mode i) 


Mode i 

Parameter Pi 

1 

0.603527579439Z? + 00 

2 

0.4230902897687? + 00 

3 

0.2323644891231? + 01 

4 

0.6896051022711? - 02 

5 

0.1758143354571? - 01 

6 

0.3391290367501? - 02 

7 

0.5186700362261? - 04 

8 

0.1768623404521? - 03 

9 

0.4327256540097? - 03 
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CHAPTER 5 


FINITE ELEMENT WING MODELING USING NASTRAN 
FOR AEROELASTIC SENSITIVITY CALCULATIONS 

5.1 OVERVIEW 

In previous chapters, to obtain sensitivity of aeroelastic response of a wing 
to shape and modal parameters, a Rayleigh-Ritz based equivalent plate model was 
used for analyzing the wing structure. Though equivalent plate representation 
in conjunction with global Ritz analysis techniques are relatively inexpensive for 
mathematical modeling of wings to study aeroelasticity, finite element methods 
[91-93] are being used in the industry for structural modeling of aircraft wings. 

Several commercial finite element codes like ABAQUS, ANSYS, EAL, IDEAS 
and NASTRAN exist for finite element analysis of structures. The current work 
focusses on modeling the wing structure using NASTRAN [94]. One of the features 
of using FEM is that any planform shape and wing cross-section could be analyzed 
as opposed to the Rayleigh-Ritz formulation based equivalent plate model used 
in Chapters 3 and 4, where only a trapezoidal skewed configuration of uniform 
thickness and layerwise construction could be analyzed. The wing skins, spars 
and ribs of an airplane wing can be modeled using shell and beam elements in 
NASTRAN. The free vibration modes of the wing can be fed into the unsteady 
aerodynamics code to generate the aerodynamic forces for aeroelastic analysis. 
The sensitivity of flutter speed of the wing to shape and modal parameters is 
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computed using a combination of the central difference scheme and the automatic 
differentiation software ADIFOR. 

5.2 NASTRAN ELEMENTS USED IN THIS STUDY 

In this study, the CQUAD4 and CTRIA3 shell elements, the CBAR beam 
elements and the CSHEAR shear panel elements are used. A description of these 
elements is given below. 

5.2.1 CQUAD4 SHELL ELEMENT 

The CQUAD4 element is one of the isoparametric quadrilateral shell elements 
in NASTRAN with optional coupling of bending and membrane stiffness. The 
element coordinate system for quadrilateral shell element is shown in Fig. 5.1. 
The integers 1, 2, 3 and 4 refer to the order of the connected grid points on the 
connection entries defining the elements. The angle, THETA, is the orientation 
angle for material properties. The CQUAD4 element is a 4-noded, 24 degrees of 
freedom element with 6 degrees of freedom per node (u, v, w, 6 X , Q y , 9 Z ). 

5.2.2 CTRIA3 SHELL ELEMENT 

The CTRIA3 element is one of the isoparametric triangular shell elements 
in NASTRAN with optional coupling of bending and membrane stiffness. The 
element coordinate system for triangular shell element is shown in Fig. 5.2. The 
integers 1, 2 and 3 refer to the order of the connected grid points on the connection 
entries defining the elements. The angle, THETA, is the orientation angle for 
material properties. The CTRIA3 element is a 3-noded, 18 degrees of freedom 
element with 6 degrees of freedom per node («, v, w, 9 X , d y , 9 Z ). 
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5.2.3 CBAR BEAM ELEMENT 


The CBAR element is a one-dimensional bending element which is pris- 
matic, and for which the elastic axis, gravity axis, and shear center all coincide. 
The bar element coordinate system is shown in Fig. 5.3. The CBAR element 
is a 2-noded, 12 degree of freedom element with 6 degrees of freedom per node 
(a, v, w, 9 X , 6 y , 9 Z ). 


5.2.4 CSHEAR SHEAR PANEL ELEMENTS 

The CSHEAR shear panel element is a two-dimensional structural element 
that resists the action of tangential forces applied to its edges, and the action 
of normal forces if effectiveness factors are used. The use of proper effectiveness 
factors to account for the normal forces on the elements are described in section 
5.3.2 where shear panels are used to model the spar and rib webs. The structural 
and non-structural mass of the shear panel is lumped at the connected grid points. 
The element coordinate system for a shear panel is shown in Fig. 5.4a. The labels 
Gl, G2, G3 and G4 refer the order of the connected grid points. The element 
forces are shown in Fig. 5.4b, which consist of the forces applied to the element 
at the corners in the direction of the sides, kick forces at the corners in a direction 
normal to the plane formed by the two adjacent edges, and shear flows (force per 
unit length) along the four edges. 
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5.3 NUMERICAL EXAMPLES 


In order to validate the finite element modeling, static analysis of a box-beam 
is performed using NASTRAN and the results compared with Euler-BernoulH solu- 
tion. Both static and dynamic analysis of a swept-back wing are carried out using 
NASTRAN and the wing deflections and frequencies are compared with previously 
published results. 

5.3.1 FINITE ELEMENT BOX-BEAM MODEL USING NASTRAN 

As a first example, a cantilevered box-beam shown in Fig. 5.5 was analyzed. 
The box-beam is modeled with CQUAD4 shell elements and CBAR beam elements. 
A total of 70 plate elements and 60 beam elements were used for the discretization. 
The dimensions of the box-beam are as follows: 

Thickness of the plate elements, t == 5 mm 

Length of the beam, L — 4m 

Width of the beam, b = 0.8 m 

Height of the beam, h = 0.4m 

Area of the beam elements, Ab = 900mm 2 

Young’s Modulus, E = 210GPa 

Poisson’s ratio, v — 0.3 

Moment of inertia, I = 586 x 10 -6 m 4 

The box-beam was subjected to two loadings: (i) a tip load of P = 300 iV, 
and (ii) a distributed load of q = 80 N/m (corresponding to a uniform pressure load 
of 100 N lm 2 on the beam 0.8m wide). The tip deflection results from NASTRAN, 
at nodes 1 to 6 (see Fig. 5.5) on the tip of the beam, were compared with Euler- 
Bernoulli beam-theory solution for a beam given below: 
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For a tip load, 


For a distributed load, 


w = 


PL 3 
3 El 


w = 


gX 4 

8EI 


(5.1) 


(5.2) 


The results are shown in Fig. 5.6. There is good agreement between the finite 
element solution and the beam-theory solution. 


5.3.2 STATIC AND DYNAMIC ANALYSIS OF A SWEPT WING 

In order to validate the modeling, a swept back wing model reported in Ref. 
95 published for AGARD is analyzed. It is an all aluminum wing, swept back 
by 30 deg. and has 5 identical spars and 3 identical ribs bonded to the top and 
bottom cover skins. The dimensions of the wing and its cross-section properties 
are given in Fig. 5.7. The finite element discretization of the wing is given in Fig. 
5.8. 

The wing is subjected to a tip load of 1 lb. at the trailing edge as shown 
in Fig. 5.9. The wing skins are modeled using CQUAD4 and CTRIA3 shell 
membrane-bending elements. The spar and rib webs are modeled using CSHEAR 
panel elements and the spar and rib caps axe modeled using CBAR beam elements. 
The computed deflections from NASTRAN of the front and rear spars along the 
span are plotted in Fig. 5.9. The measured deflection from Ref. 95 is also shown. 
The results agree fairly well. 

Next, the natural frequencies of the wing are calculated using NASTRAN. 
The wing skins are modeled using CQUAD4 and CTRIA3 shell membrane-bending 
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elements and the spar and rib caps are modeled using CBAR beam elements as 
before. However, the spar and rib webs are modeled in two different ways: (i) 
using CSHEAR shear panel elements and (ii) CQUAD4 shell membrane-bending 
elements. The natural frequencies of the wing are compared in Table 5.1 to the 
results reported in Ref. 96 which were computed using ELFINI finite element 
package. It should be noted that in the finite element model reported in Ref. 96, 
the webs act only in shear. The results obtained using shear panels for the webs 
compare well with ELFINI results. When shear panels are used to model the webs, 
it is assumed that the webs act only in shear and the required axial stiffness at 
the edges are provided to support the corner forces. However, these webs also 
contribute to the stiffness of the structure in a small way. Hence, the spar and 
rib webs are modeled using CQUAD4 shell membrane-bending elements and the 
first five frequencies are given in Table 5.1. The CQUAD4 element for the webs 
gives stiffer solutions as can be seen from Table 5.1 (the frequencies of transverse 
vibration are about 2-3 % higher, especially in the bending modes, and the inplane 
frequency is also higher). 

When using the CSHEAR element in NASTRAN (see Fig. 5.4), appropriate 
factors Fi and . 2*2 have to be specified for providing extensional stiffness along 
the sides of the element to handle the corner forces which are parallel to the 
sides of the element. Factor Fi is the effectiveness factor for extensional stiffness 
along edges 1-2. and 3-4 and 2*2 is the effectiveness factor for extensional stiffness 
along edges 2-3 and 1-4. The effective extensional area is defined* by means of 
extensional rods on the perimeter of the element. Thus, if F\ = 1.0, the panel 
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is fully effective for extension in the 1-2 and 3-4 direction. The areas of the rods 
on the edges 1-2 and 3-4 are set equal to (0.5 F\ T w\), where w\ is the average 
width of the panel and T is the thickness of the panel. Proper selection of the 
effectiveness factors is important when using the shear panel elements. Three 
different cases were examined using the effectiveness factors. If the default values 
of F\ = i *2 = 0.0 are used, then the computed frequencies axe very small, since 
the shear panel elements have no axial stiffness along the edges. If the panels are 
made fully effective for extension in both directions, then the solution obtained is 
much stiffer than that obtained using CQUAD4 shell membrane-bending elements 
for the webs. The proper choice is to make the panel (for spar and rib webs) fully 
effective for extension in the transverse direction of the wing and let the spar and 
rib caps carry the load in the other direction. This gives frequencies which agree 
well with the frequencies reported in Ref. 96 using ELFINI. 

5.4 SENSITIVITY OF FLUTTER SPEED 

The structural modeling of the wing is done using NASTRAN. The airfoil 
shape is generated by transforming a circle into an airfoil using the Joukowski 
transformation, so that shapes of varying thickness and camber can be obtained 
by varying the parameters. Finite element discretization of the wing skins, spars 
and ribs is done and a free vibration analysis of the wing is performed. The 
natural frequencies and mode shapes obtained are used to generate the generalized 
aerodynamic forces required for flutter analysis using FAST. The sensitivity of 
flutter speed to shape and modal parameters is computed using a combination of 
central difference scheme and ADIFOR. 
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5.4.1 AIRFOIL COORDINATES FOR THE WING 


It is beneficial if the airfoil shape generation for the wing cross-section is 
parameterized so that shapes of varying thickness and camber can be generated. 
In order to generate the co-ordinates of the wing cross-section, a transformation 
of a circle into an airfoil using the Joukowslci transformation is carried out. The 
equation of the circle in the (-plane shown in Fig. 5.10 is 

(Z-P ) 2 +(v~<l ) 2 =« 2 (5-3) 


with center at p — p + iq and radius a. 
For points on the circle, 


f == p + acosd 
r? = q + asind 


0 < d < 2tt 


Therefore, 


( = £ 4. irj = (p 4- acosd ) + i{q + asind ) 


(5.4) 


(5.5) 


The circle is transformed into an airfoil in the 2 -plane using the Joukowski 
transformation 

* = C + ~r (5.6) 

Substituting for ( from equation(5.5) and denoting 


D = p 2 + q 2 + a 2 + 2 apcosd 4- 2 aqsind (5-7) 

we have 

-=(f + -cosd)( 1 4- — ) + 1 {- C + -Sind) (l - — ) (5.8) 

87 



A number of different airfoil shapes can be generated by varying the values 
of p and q. Choosing p = —c/12 and q = 0, the airfoil generated is shown in Fig. 
5.10. 

5.4.2 FREE VIBRATION MODES OF THE WING 

Before performing a flutter analysis using the lifting-surface unsteady aero- 
dynamic program FAST, the free vibration analysis of the wing structure has to 
be performed. The finite element discretization of the wing is shown in Fig. 5.11, 
with the coordinates of the upper and lower skins of the wing obtained from equa- 
tion (5.8). The airfoil thickness at any chord location varies linearly from the root 
of the wing to the tip. The wing is modeled using 4 spars (placed at 0.05c, 0.2c, 
0.5c and 0.8c location) and 10 ribs (placed equidistant along the span). A total 
of 190 CQUAD4 and CTRIA3 shell elements and 200 CBAR beam elements were 
used to model the wing structure. The dimensions of the wing are as follows: 

Wing span, L — 5 m 

Area, .4 = 7.5 m 2 

Root chord, C r — 2 m 

Tip chord, Ct — 1m 

Wing skin thickness, t = 3 mm 

Area of beam elements, Aj = 75mm 2 

Young’s Modulus, E = lOGPa 

Poisson’s ratio, v = 0.3 

Material density, p = 2700 kg/m 3 

The natural frequencies and the mode shapes of the wing are obtained from 
NASTRAN and the first six frequencies and the corresponding free vibration modes 
of the wing are given in Fig. 5.12. 
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5.4.3 SENSITIVITY TO SHAPE AND MODAL PARAMETERS 


In order to perform the aeroelastic analysis using FAST, six vibration modes 
of the wing are used. The modal information from NASTRAN finite element 
analysis is input into FAST to generate the flutter speed. The sensitivity analysis 
is carried out in two steps — a finite difference calculation followed by automatic 
differentiation. Since the new coordinates of the finite element nodal points are 
generated by perturbing the value of the shape parameter and these are input 
into NASTRAN in an 8-column format, it was observed that the precision of the 
numbers was not good enough which could lead to errors during finite differencing. 
So the central difference method was preferred over the forward finite difference 
method. One of the shape parameters of the wing is peturbed by a small positive 
value and the free vibration modes of the new wing configuration are obtained from 
NASTRAN. Again, the shape parameter is perturbed by a small negative value 
from its nominal value and modal analysis is carried out. It should be noted that 
the repeated finite element calculation is not very expensive, because the analysis 
can be performed to yield only a certain number of modes which are required for 
the flutter analysis, rather than compute all the frequencies. The derivatives of the 
modal frequency and mode shapes are calculated using a central difference scheme 

df f{x + Ax) - f{x - Ax) . . 

dx 2Ax (5 ' 9) 

When the modal information is read into FAST for flutter analysis, the corre- 
sponding derivative information is also read in. These derivative values are then 
propagated and the derivatives of the subsequent active variables are calculated 
using automatic differentiation. 
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The sensitivity of flutter speed to shape parameters is plotted in Figs. 5.13- 
5.16 for M = 0.6 and p az r — 0.8 kg/m 3 along with the flutter speed calculated 
from a reanalysis by changing one shape parameter at a time. The shape sensitiv- 
ity derivatives give the trend that can be expected in the flutter speed variation for 
small changes in shape parameters about the baseline configuration. The deriva- 
tives of the flutter speed with respect to modal parameters are also calculated using 
ADIFOR, and the results are given in Figs. 5.17-5.19. It can be seen that the 
higher modes do not actively participate in flutter as evidenced by the sensitivity 
of the flutter speed to the modal parameters of the higher modes. 
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Fig. 5.1 CQUAD4 shell element 



Fig. 5.2 CTRIA3 shell element 
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Fig. 5.6 Static deflection of the box-beam 
for two load cases 
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Fig. 5.7 Planform dimensions and cross-section properties of the wing 
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Fig. 5.8 Finite element discretization of the AGARD swept back wing model 
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Fig. 5.9 Deflection of front and rear spar for 
a tip load of 1 lb. at the rear spar. 
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Fig. 5.13 Flutter speed vs Aspect ratio (M=0.6) 

(AR=3.33, Area=7.5nrf, TR=0.5, Sweep=15°) 
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Fig. 5.16 Flutter speed vs Sweep angle (M=0.6) 

(AR=3.33, Area=7.5m 2 , TR=0.5, Sweep=15°) 






Fig. 5.17 Flutter speed sensitivity to natural frequency (M=0.6) 

in (m/s)/(rad/s) 

(AR=3.33, Area=7.5m 2 , TR=0.5, Sweep=1& } ) 
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Fig. 5.18 Flutter speed sensitivity to generalized mass (M=0.6) 

in (m/s)/(kg) 

(AR=3.33, Area-7.5rrf, TR=0.5, Sweep^S 3 ) 



( 1 , 1 ) ( 2 , 1 ) ( 3 , 1 ) ( 6 , 1 ) ( 3 , 2 ) ( 3 , 3 ) ( 1 , 4 ) ( 6 , 6 ) 

GAF 

Fig. 5.19 Flutter speed sensitivity to nondimensional generalized 
aerodynamic forces (GAF) (M=0.6) in (m/s) 

(AR=3.33, Area=7.5m 2 , TR=0.5, Sweep=1&) 
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Table 5.1 Natural frequencies of the AGARD Swept Wing Model 


Mode 

Number 

Frequencies (Hz.) 

Present (CSHEAR) 

Livne[96] 

Present (CQUAD4) 

1 (1st Bending) 

116.8 

115.6 

120.8 

2 (Inplane) 

318.0 

317.6 

353.0 

3 (1st Torsion) 

415.4 

418.4 

416.9 

4 (2nd Bending) 

579.5 

576.4 

593.1 

5 (2nd Torsion) 

1082.3 

1086.0 

1094.3 
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CHAPTER 6 


SHAPE SENSITIVITY FOR OPTIMIZATION 
WITH AEROELASTIC CONSTRAINTS 

6.1 OVERVIEW 

An equivalent plate model which incorporates first-order shear deformation 
theory is examined so it can be used to model thick wings, where shear deforma- 
tions are important. The natural frequencies of a thick laminated composite plate 
obtained using first-order shear deformation theory are compared to those obtained 
using classical plate theory. The sensitivity of natural frequencies to changes in 
shape parameters is obtained using ADIFOR. 

The use of shape sensitivity derivatives in a gradient-based design optimiza- 
tion subjected to aeroelastic constraints is demonstrated. A simple optimization 
effort is made towards obtaining a minimum weight design of the wing, subject to 
flutter constraints, lift requirement constraints for level flight and side constraints 
on the planform parameters of the wing. The IMSL subroutine NCONG, which is 
an implementation of successive quadratic programming is used for this purpose. 

6.2 STRUCTURAL MODEL 

In the previous chapters, equivalent plate models based on classical plate 
theory and finite element models were used for structural modeling. The first- 
order shear deformation theory which incorporates transverse shear developed by 
Kapania and Lovejoy [97] is used for structural modeling in this chapter. 
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The plate displacements u, v and to are given by 

u = u°{x,y,t) + Z(j>x(x,y,t) 

v = v°(x,y,t) + z<j> y (x,y,t ) ( 6 . 1 ) 

w = w°(x,y,t ) 

where u°, v° and w° are midplane displacements and <j> x and 4> y are rotations about 
the y and x axes, respectively. 

The original (x, y) coordinate system is transformed to the ( 77 , £) coordinate 
system, and a Rayleigh- Ritz formulation using Chebyshev polynomials to represent 
the plate displacements is adopted. The mid-plane displacements and rotations 
are represented in a series as 

I J 

“” = E 

i=Q j—Q 
K L 

b “ = EE s u (mMT,(o 

k—0 1—0 
M N 

w ° = E E p mn(t)T m (rj)T n (0 (6.2) 

m=0 n=0 

P Q 

p = 0 g= 0 

R S 
r=0 5=0 

and the Chebyshev polynomials T; are given in equations(3.1). 

It has been shown by Kapania and Lovejoy [97] that for free vibration of the 
plate, the eigenvalue problem can be obtained as 

[K - XM]{x} = 0 (6.3) 
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where [A”] and [M] are the stiffness and mass matrices and the eigenvector {or} 
contains coefficients of the polynomial representing the displacement functions, 
i.e., 

} ~ (Roo > Roi 7 • ■•Rij] S 00 7 •••Ski] P DO 7 ■••Pmn ] V0O7 V pq ] ^ 00 7 ( 6 . 4 ) 

The plate boundary conditions are handled by using springs of appropriate 
magnitude at the boundaries. For modeling cantilever wings, linear and rotational 
springs of large magnitude are placed at the root to satisfy the clamped wing 
boundary condition. Details of the formulation are given in [97]. 

Both inplane modes of vibration and rotary inertia terms are included in [97]. 
However, in this study, considering only transverse vibrations and neglecting the 
rotary inertia, the coefficients corresponding to the u°, v°, (j> x and <p y displacements 
are condensed out using a static condensation. 

6.3 EFFECT OF TRANSVERSE SHEAR ON VIBRATION AND 
AEROELASTIC CHARACTERISTICS 

The natural frequencies of a laminated composite plate clamped at one end 
and free at the other three edges, are computed using this structural model in- 
corporating transverse shear and compared with the classical plate model used 
in the earlier studies. The sensitivity of natural frequencies to changes in shape 
parameters of the plate is studied. The influence of transverse shear on the flutter 
characterisitics of a wing is also examined. 
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6.3.1 EFFECT OF TRANSVERSE SHEAR ON NATURAL FREQUENCY 


In order to get an estimate of the effect of transverse shear on the free vi- 
bration characteristics, the plate vibrations are analyzed using both the classical 
plate theory and the first-order shear deformation theory. The results for sym- 
metrically laminated, skew, cantilever plates made of Glass/Epoxy are presented 
in Tables 6.1 and 6.2. The properties for Glass/Epoxy are E\ = 38.61 GPa , 
E 2 — Ez = 8.27 GPa , G \2 — G 13 = 4.14 GPa, G 23 — 3.35 GPa, v \2 — v\z — 0.26, 
1/23 = 0.49 and p = 2546.54 kg/m 3 . The results axe given for a plate with a 
sweep angle of A = 15°, aspect ratio of 3.111, taper ratio of 0.5 and area of 
63 in 2 (406.45 cm 2 ). The thickness of the plate is varied from 0.14 in. (3.556 mm.) 
for the thin [302 /0] s laminate to 0.84 in. (21.336 mm.) for the thick [30s /0] 6 3 lami- 
nate. It is seen from Tables 6.1 and 6.2 that the frequencies predicted by including 
shear deformations are lower than those predicted by classical plate theory. As the 
thickness of the plate is increased, notable differences can be seen in frequencies 
obtained by the two theories. 

6.3.2 SHAPE SENSITIVITY RESULTS 

The sensitivity of the natural frequencies of the plate to shape parameters, 
namely aspect ratio, area, taper ratio and sweep angle was computed using au- 
tomatic differentiation (AD IFOR). The [ 3 O 2 /0] io.s laminated plate with 60 plies, 
made of Glass/Epoxy is used for the study. The sensitivity of the first three nat- 
ural frequencies of the plate with respect to shape parameters is shown in Figs. 
6. 1-6.4. The free vibration frequencies obtained by perturbing one shape parame- 
ter at a time from the baseline configuration are shown by the solid line curve. The 
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dotted line gives the sensitivity prediction at the baseline configuration. The sen- 
sitivity derivative obtained forms a tangent to the frequency curve at the baseline 
configuration. 

It is seen that there is a drop in the vibration frequencies of the first three 
modes as the aspect ratio, area and taper ratio of the plate are increased. The 
first two natural frequencies decrease with increasing sweep angle as seen from 
Fig. 6.4, but the third natural frequency first increases with sweep angle and then 
decreases, about the baseline value of 15° sweep. By performing one sensitivity 
calculation at the baseline analytically, this method gives a linear approximation to 
the frequency curve for changes in the wing shape parameters about the baseline. 
These shape sensitivity derivatives could be very useful for a gradient-based design 
optimization of thick plates with frequency constraints. 

6.3.3 EFFECT OF TRANSVERSE SHEAR ON FLUTTER SPEED 

Studies on the effect of transverse shear on flutter speed of a wing is of 
great interest, especially when dealing with thick composite wings for which shear 
deformation effects are not negligible (see Karpouzian and Librescu [98]). In this 
study, the first-order shear deformation theory (FSDT) is used for the structural 
model and a lifting-surface unsteady aerodynamics (FAST) is used for generating 
the generalized aerodynamic forces required for the flutter calculation. The flutter 
characteristics of the wing shown in Figure 3.2 which has previously been analyzed 
in Chapter 4 using a classical laminated plate theory (CLPT), is studied here by 
including the effects of transverse shear. Figures 6.5-6. 8 show the variation of 
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the flutter speed with respect to aspect ratio, area, taper ratio and sweep angle 
of the wing, as predicted by CLPT and FSDT. It is evident that the inclusion of 
transverse shear in the structural model has an effect of reducing the flutter speed 
of the wing. 

6.4 OPTIMIZATION WITH AEROELASTIC CONSTRAINTS 

Multidisciplinary optimization is of great interest in the aerospace industry 
and efforts are being made towards integrated wing synthesis to produce better 
designs. The designs have to satisfy a number of constraints of which flutter con- 
straints form an important subset. The changes in shape of the wing affect both its 
structural and aerodynamic response. This is an effort to obtain a minimum weight 
design of the wing, subject to flutter constraints, lift requirement constraints and 
constraints on the planform parameters of the wing. It should be noted that more 
constraints can be added to this optimization problem if required. The gradients 
required for the optimization are computed through a sensitivity analysis of the 
aeroelastic response of the wing using automatic differentiation (using ADIFOE). 

6.4.1 FORMULATION OF THE OPTIMIZATION PROBLEM 

The optimization problem can be stated as follows: 

Minimize 

= pmStg 

subject to the constraints 

V flutter ^ 1-2 Vf light 
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I = \pVj, igh ,SC L = 

AR > 5.0 
0.5 < TR < 1.0 

0° < A < 30° (6.5) 

where W w is the weight of the wing, p m is the density of the material, S is the 
area of the wing, t is the thickness, g is the acceleration due to gravity, V is the 
velocity, L is the lift, p is the density of air, Cl is the lift coefficient, W ac is the 
weight of the aircraft, AR is the aspect ratio, TR is the taper ratio and A is the 
quarter-chord sweep angle. 

The wing is designed so that the flutter speed is atleast 20% more than the 
maximum flight speed. Also the wing is capable of producing sufficient lift required 
for cruise. The constraint on the aspect ratio ensures that wings with short spans 
and long chords are not an outcome of the optimization. The constraints on taper 
ratio and sweep angle ensure that wings that conform to practical designs are 
generated. 

6.4.2 METHOD OF SOLUTION 

There are several methods to perform this nonlinear constrained optimization 
[99] . The gradient projection method is based on projecting the search direction 
into the subspace tangent to the active constraints. In the method of feasible 
directions, the concept is to stay in the feasible domain, move in a direction which 
reduces the objective function and stay away from the constraint boundaries. The 
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program by Vanderplaats, CONMIN [100] is an implementation of the method of 
feasible directions. Another method which uses successive quadratic programming 
is implemented in the IMSL subroutine, NCONG [101]. The method, based on the 
iterative formulation and solution of quadratic programming subproblems, uses a 
quadratic approximation of the Lagrangian and linearization of the constraints. 
Though the solution of the quadratic programming direction seeking problem is 
cumbersome, it often leads to faster convergence. 

The optimization problem can be written as 


( 6 . 6 ) 


Min /(x) 

subject to gj(x) = 0, for j = 1, . . . , m e 

gj(x) > 0, for j = m e + 1, . . , , m 

We seek the direction d as the solution of the following quadratic program- 
ming problem: 


su 


Min id T B*d + V/(Xfc) T d 
bject to Vgj(x k ) T d + gj(x k ) = 0, j = 1, . . . , 


m e 


(6.7) 


v J>(x/i) :r <i + s,(x t ) >0, j = m, + 
where Bfc is a positive definite approximation of the Hessian, and x k is the current 

iterate. If dfc is the solution to the subproblem, a line search is used to find the 

new point x k +i, 

Xfc+i = x fc + Ad*, A € (0, 1] (6.8) 


Bfc is updated according to the modified BFGS formula [102]. It should be noted 
that this algorithm can have infeasible points during the solution process. 
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6.4.3 OPTIMIZATION RESULTS 


The optimization problem described in equations(6.5) was solved using the 
IMSL routine NCONG. The wing is shown in Fig. 3.2. The wing skins axe made 
of 0 deg. laminated Graphite/Epoxy (T300/iVo208) with the following material 
properties: E\ — 181 x 10 9 Pa, E z = 10.3 x 10 9 Pa, E z — 10.3 x 10 9 Pa, 
ut 2 = 0.28, i/ia = 0.28, u 2Z = 0.49, G n = 7.17 x 10 9 Pa, G 1Z = 7.17 x 10 9 Pa, 
G23 = 5.82 x 10 9 Pa and p m = 1600 kg/m 3 . The analysis is carried out at a Mach 
number of 0.6. The problem was further simplified by fixing the taper ratio and 
sweep angle of the wing at 0.5 and 30°, respectively, so the design variables are 
the aspect ratio and the area of the wing. The following values were chosen: 


Vf light — 100 m/ S 


p = 0.8 kg/m 3 
C L = 0.7 


(6.9) 


W ac = 100000 N 

The gradients of the aeroelastic constraints are obtained by performing sensi- 
tivity analysis of the aeroelastic response with respect to the planform parameters 
using the automatic differentiation package, ADIFOR. The results of the optimiza- 
tion for each iteration are given in Table 6.3. Since the material and thickness of 
the wing are not changing, the wing weight is proportional to the area of the wing. 
It can be seen that at the 3rd iteration, AR — 5.3913 and Area — 17.857 m 2 , 
corresponding to Vfi u u e r — 134.19 m/s. This value of the wing area satisfies the 
equality constraint, which means it is not possible to obtain a reduction in the 
value of the objective function by reducing the area since the equality constraint 
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will then be violated. The minimum weight of the wing, subject to the applied 
constraints is W w = [17.857 p m tg] N. 

In this optimization, the equality constraint appears to be the driving con- 
straint and the inequality constraint on the flutter speed is not active. Even though 
both the constraints are satisfied and the minimum weight as desired by the ob- 
jective function is achieved, there is a big margin on the flutter speed constraint 
(134.19 — 120.00 = 14.19 m/s), which suggests that we could have other designs 
with higher aspect ratio which will have the same optimum weight value and make 
both the constraints active. Since the objective function is only a function of area, 
once its minimum is attained by the satisfaction of the equality constraint, the 
optimization algorithm does not find a direction in which further minimization of 
the objective function can be accomplished. This either results in the inequality 
constraint being satisfied with a wider margin as in this case, or in the case of 
a violated constraint, it becomes necessary to optimize on the single design vari- 
able namely the aspect ratio to meet the inequality constraint. This is avoided 
by defining a composite objective function by adding to the area of the wing the 
reciprocal of the difference between the aspect ratio and its lower bound, weighed 
with a parameter p as shown. 


/ 


S + p 


( AR - ARib) 


( 6 . 10 ) 


where / is the new objective function, S is the wing area, AR is the aspect ratio 
and ARib is its lower bound. With this newly defined objective function, the 
results from the optimization are given in Table 6.4 for a value of p = 10.0. It 
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can be seen from Table 6.4 that from the 3rd iteration on, the Area = 17.857 m 2 , 
which is driven by the equality constraint. However, the optimization algorithm is 
successful in its search for a direction which minimizes the objective function and 
proceeds until both the constraints are active, giving a final value of AR = 6.5994 
and Area — 17.857 m 2 with Vfi utter = 120.00 m/s. The minimum weight of the 
wing is W w = [17.857p m t#] N. 

This is a simple two constraint optimization problem with two design vari- 
ables. It should be noted that the material properties and skin thicknesses of the 
wing could also be used as design variables for a full-scale optimization with aeroe- 
lastic constraints; however a much simpler problem was chosen here to demonstrate 
the use of shape sensitivity derivatives in optimization. 
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Fig. 6.1 First three natural frequencies vs Aspect ratio 
for a thick plate 

(AR=3.1 1 1, Area=0. 040645m 2 , TR=0.5, Sweep=15°) 
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Fig. 6.2 First three natural frequencies vs Area 
for a thick plate 

(AR=3. Ill, Area=0.040645rrf, TR=0.5, Sweep=15°) 
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Fig. 6.3 First three natural frequencies vs Taper ratio 
for a thick plate 

(AR-3.111, Area=0.040645rrf, TR-0.5, Sweep=15°) 
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Fig. 6.4 First three natural frequencies vs Sweep angle 

for a thick plate 

(AR=3.111, Area=Q, 040645 rrf, TR=0.5, Sweep=15°) 
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Fig. 6.5 Effect of transverse shear on fiutter 
[Flutter speed vs Aspect ratio (M=0.6)] 

(Area=20m 2 , TR=0.5, Sweep=30 p ) 
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Fig. 6.6 Effect of transverse shear on flutter 
[Flutter speed Vs Area (M=0.6)] 

(AR=10, TR=0.5, Sweep=3Cf>) 
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Fig. 6.7 Effect of transverse shear on flutter 
[Flutter speed vs Taper ratio (M=0.6)] 

(AR=10, Area=20nrf, Sweep=30 p ) 
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Fig. 6.8 Effect of transverse shear on flutter 
[Flutter speed vs Sweep angle (M=0.6)] 

(AR=10, Area=20rrf, TR=0.5) 




Table 6.1 Comparison of natural frequencies predicted by the two theories 
for the 6-ply, 12-ply and 18-ply laminates 


Mode 

Number 

Natural frequencies (Hz.) 

[30 2 /0] a 

(6-ply) 

[30 2 /0] 23 (12-ply) 

[30 2 /0] 

3 * (18-ply) 

CLPT 0 

FSDT 6 

CLPT 

FSDT 

CLPT 

FSDT 

1 

16.16 

16.14 

34.49 

34.42 

52.76 

52.57 

2 

80.19 

79.95 

170.23 

168.58 

259.90 

255.22 

3 

118.30 

117.30 

235.47 

230.45 

351.96 

339.61 

4 

209.92 

208.80 

445.04 

437.82 

679.59 

658.99 

5 

316.95 

313.86 

636.14 

619.30 

952.52 

908.13 


“ CLPT (Classical laminated plate theory) 

6 FSDT (First-order shear deformation theory) 


Table 6.2 Comparison of natural frequencies predicted by the two theories 

for the 24-ply and 36-ply laminates 


Mode 

Number 

Natural frequencies (Hz.) 

[30 2 /0] 43 (24-ply) 

[30 2 /0] 63 (36-ply) 

CLPT 

FSDT 

CLPT 

FSDT 

1 

71.08 

70.61 

107.57 

106.30 

2 

349.77 

339.71 

528.86 

501.06 

3 

468.91 

444.58 

701.67 

640.29 

4 

914.84 

870.50 

1383.71 

1259.17 

5 

1270.76 

1178.56 

1902.68 

1658.03 
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Table 6.3 Optimization results for the minimum weight wing design 


Iteration 

Aspect ratio 

Area (m 2 ) 

Flutter speed (m/s) 

1 

10.000 

20.000 

94.30 

2 

5.3913 

18.551 

130.90 

3 

5.3913 

17.857 

134.19 


Table 6.4 Optimization results for the minimum weight wing design 

(Composite objective function) 


Iteration 

Aspect ratio 

Area (m 2 ) 

Flutter speed (m/s) 

1 

10.000 

20.000 

94.30 

2 

7.7882 

19.304 

105.56 

3 

6.1784 

17.857 

124.21 

4 

6.5694 

-17.857 

120.28 

5 

6.5994 

17.857 

120.00 
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CHAPTER 7 

CONCLUDING REMARKS 


This study has examined the flutter characteristic of a typical section in 
subsonic compressible flow. Indicial response functions are used for the normal 
force and pitching moment coefficients of the two degree of freedom airfoil and 
the unsteady aerodynamics is represented by a state-space model. The aeroelastic 
equations are solved as an eigenvalue problem to determine the onset of aeroelastic 
instability. The sensitivity of the flutter speed of the typical section with respect 
to its mass and stiffness parameters, namely, mass ratio, static unbalance, radius 
of gyration, bending frequency and torsional frequency are calculated analytically. 
The sensitivity derivatives show good agreement with finite difference derivatives. 
The aeroelastic equations are also integrated with respect to time at different 
values of freestream speed, to observe the aeroelastic phenomena in real time. The 
stable, neutrally stable and unstable nature of the aeroelastic response can be 
seen at speeds below, at and above the flutter speeds, respectively, for the typical 
section. 

The above aerodynamics for a typical section was extended to develop a strip- 
theory formulation representing the unsteady aerodynamic forces on a wing. The 
structural modeling is done using classical plate theory and is based on a Rayleigh- 
Ritz formulation using Chebyshev polynomials for the wing displacements. The 
structural equations are coupled with the time-domain aerodynamic equations to 
formulate the aeroelastic equations as a set of first-order ordinary differential equa- 
tions. These equations are solved as an eigenvalue problem to determine the critical 
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speed of the wing. The natural frequencies and flutter speeds are compared with 
previously published experimental values obtained from wind tunnel tests and the 
results agree fairly well. The sensitivity of divergence and flutter speeds to shape 
parameters, namely, aspect ratio, area, taper ratio and sweep angle are computed 
analytically. These derivatives have been accurately evaluated and they form a 
tangent to the critical speed curve at the baseline configuration. These shape sen- 
sitivity derivatives give a linear approximation to the critical speed curve about 
the baseline configuration and will be useful for preliminary design or an optimiza- 
tion with aeroelastic constraints. The time integration of the aeroelastic equations 
shows the real-time aeroelastic response (tip deflection as a function of time) of the 
wing when operating at speeds at the verge of instability. Particularly notable is 
the response of the wing configuration where the divergence speed and the flutter 
speed are close to each other. 

Flutter analysis of the wing is also carried out using a lifting-surface subsonic 
kernel function aerodynamic theory (FAST) and an equivalent plate structural 
model. The generalized aerodynamic forces are obtained for a fixed number of 
free vibraton modes of the wing, for the specified Mach number and spanning the 
range of reduced frequencies. The flutter speed is obtained using a V-g type of 
solution. The novel method of automatic differentiation using AD IFOR was suc- 
cessfully implemented to generate exact derivatives of the flutter speed as obtained 
from the discrete system, with respect to shape parameters of the wing. A good 
sensitivity prediction is obtained using ADIFOR. Also, based on the sensitivity 
of flutter speed to modal parameters, namely, natural frequency, generalized mass 
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and generalized aerodynamic forces, computed using ADIFOR, one could make 
a judicious choice of the number of modes to be used in the aeroelastic analysis 
for reasonably estimating the flutter speed. These derivatives give a qualitative 
estimate of the modes that actively participate in flutter. In order to obtain a 
quantitative estimate of the contribution of a particular mode to flutter, a pa- 
rameter was formed by constructing a logarithmic derivative of the sensitivity to 
modal parameters and summing up the absolute values of these derivatives for each 
mode. The finitely vanishing value of the' parameter for higher modes indicates 
that even if these modes are not used for the flutter calculation, the results would 
be affected only to an insignificant amount. 

Finite element modeling of the wing is done using NASTRAN so that wing 
structures made of spars and ribs and top and bottom wing skins could be analyzed. 
The modeling is validated for static and dynamic analysis using a box-beam and 
an AGARD swept-back wing model. A good agreement with previously published 
results is obtained. The airfoil shape for the wing cross-section was generated by 
transformation from a circle using the Joukowski transformation, so it could be 
parameterized. The free vibration modes of the wing obtained from NASTRAN are 
input into FAST to compute the flutter speed. The derivatives of flutter speed with 
respect to shape parameters are computed using a combination of central difference 
scheme and ADIFOR and the sensitivity to modal parameters are calculated using 
ADIFOR. A fairly good sensitivity prediction is obtained. It is seen that the higher 
vibration modes of the wing do not actively participate in flutter. 

An equivalent plate model which incorporates first-order shear deformation 
theory is then examined so it can be used to model thick wings, where shear de- 
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formations are important. The sensitivity of natural frequencies to changes in 
shape parameters is obtained using ADIFOR which would be useful in optimiza- 
tion with frequency constraints. The natural frequencies and the flutter speeds 
calculated using the first-order shear deformation theory are compared with those 
obtained using a classical laminated plate theory. It is seen that the frequencies 
and the flutter speeds drop as the transverse shear effects come into play. The 
shape sensitivity derivatives of the flutter speed are indispensable in a gradient- 
based optimization with aeroelastic constraints. A simple optimization effort is 
made towards obtaining a minimum weight design of the wing, subject to flutter 
constraints, lift requirement constraints for level flight and side constraints on the 
planform parameters of the wing. The nonlinear constrained optimization problem 
is solved using the IMSL subroutine NCONG, which uses successive quadratic pro- 
gramming. An optimum design which satisfies all the constraints is obtained. It 
should be noted that more constraints can be added to this optimization problem 
as desired. 

Further research in this area could be to incorporate a control model so it can 
be coupled with the structural and aerodynamic model already in the time-domain 
to perform aeroservoelastic studies. Supersonic kernel function aerodynamics could 
be implemented to perform sensitivity analysis in the supersonic regime. So far, 
wings with a general trapezoidal skewed planform have been analyzed. However, 
some wings have discontinuities along the leading and trailing edges and can be 
considered to be made up of multiple trapezoidal segments. Shape sensitivities of 
flutter speed of wings made up of multiple trapezoidal segments, a good example 
of which is the High Speed Civil Transport (HSCT) wing, can be analyzed. 
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APPENDIX A 


The ciif s are given by 


where 


On = — (- — ) 0 1 bi 

c 

«22 = — (~ — )^ 2 ^2 
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«33 = 
«44 = 


K a Ti 

1 


KqTl 

055 = -(hKvMTi)- 1 


06 6 = -(64 K a MTl ) 

077 = — 65/3 ( — -) 
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088 


KqmTi 


-1 


Ka(M) = [(1 - M) + tt/9M 2 (Ai6i + A 2 6 2 )] 

Kq(M) = [(1 - M) + 2TT/3M 2 (A 1 b 1 + A 2 b 2 )\ 

T - , , fA I" A3 64 + A463 

K qM (M)- [ hb ^ l _ M y 

K qM (M ) = f 7 


[l5(l-M) + 37r0M 2 & 5 J 



The Cij' 1 s are given by 


cii = C N ^)l3 2 A 1 b l 
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Cl3 




) 


M x K a Ti 
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12M'K qM Tj 

The constants are given by A\ — 0.3, A 2 = 0.7, A 3 = 1.5, A 4 = —0.5, bi = 0.14, b 2 
0.53 , 63 = 0.25 , 64 = 0.1 , 65 = 0.5. 
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Analytical derivatives 


Aspect ratio ( AR ), Area (5), Taper ratio (tr), Sweep (/i) 

The wing coordinates are (x x ,y x ), (a; 2 ,y 2 ), (£ 3 , 2 / 3 ) and ( 0 : 4 , y 4 >. 

span — vARS 

2S 

QJ* ~ . 

span( 1 -f tr) 
ct = tr cr 

x x = 0.75 cr 

x 2 = span tanA + 0.75 ct 
xz = span tanA — 0.25 ct 
£4 = —0.25 cr 

2/1 = 2/4 = 0.0y 2 = ijz = spanp = (zi -f z 2 + x z + £ 4 ) 

pp = (£2 + £ 3 ) - (£1 + £ 4 ) 
crct = cr + ct 
rt — cr — ct 

For any point ip on the quarter chord line (—1 <ip < 1) 

y = 0.5 span (1 + ip) 
x = y tanA 



span 


__ ii pp - 4£ + p) 

(£ rt — crct) 

c ~ [cr + (ct - cr)(l + ip) 0.5] cos/1 
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Local chord fc) 


dc 
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dc 
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dc 
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2 = [xj + xj + xz + x±l 


5°' 5 [1 — (1 — tr)( 1 + ip) 0.5] cos A 
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rt = cr — ct 


x=_v tanA 


drt 

Jar 

drt 

Js 

drt 

dtr 

drt 

JI 


S 0 - 5 (1 - tr) 
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0.5 rt 
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dx 
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For any parameter v, 


dp = U rt - crct)(£^f - 4 jf + jg) - (£ pp - 4s + p)(g 
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IT- = 0.0 
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dcrct \ 
dv l 
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APPENDIX B 


Main program: main.f 

program main 
c 

real x(3),y 

x(l)=1.0 

x(2)=2.0 

x(3)=3.0 

call func(x,y) 

write(*,*)x( l),x(2),x(3),y 

end 

Subroutine: func.f 

subroutine func(x,y) 
real x(3),y 

y~x(l)+2.*x(2)**2+3.*x(3)**3 

return 

end 

ADIFOR generated code: g_func.f 


subroutine g_func(g_p_, x, g_x, ldg_x, y, g_y, ldg_y) 

real x(3), y 

integer g_pmax_ 

parameter (g_pmax_ = 3) 

integer g_i_, g_p_, ldg_y, ldg_x 

real r2_p, rl_p, r9_b, r5_b, r3_v, r7_v, g_y(Idg_y), g_x(ldg_x, 
*3) 

if (g_p_ .gt. g_pmax_) then 
print *, 'Parameter g_p_ is greater than g_pmax_' 
stop 
endif 

r3_V :as x(2) * x(2) 
r2_p = 2.0e0 * x(2) 
r7_v = x(3) **(3 -2) 
r7_v = r7_v * x(3) 
rl_p = 3 * r7_v 
r7_v = r7_v * x(3) 
r5_b = 3, * rl_p 
r9_b = 2. * r2_p 
do g_i_= 1, g_p_ 

g_y(g_L) = g_x(g_i_, 1) + r9_b * g_x(g_i_, 2) + r5_b * g_x(g_ 
*i_, 3) 
enddo 

y = x( 1) + 2. * r3_v + 3. * r7_v 


return 

end 
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